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Introduction 

The  “Time-resolved  Spectral  Optical  Breast  Tomography”  research  project  aims  to  develop  a  near  real¬ 
time  three-dimensional  (3D)  spectral  tomographic  imag  ing  algorithm  with  use  of  a  transport  model  based 
on  radiative  transfer  and  light  of  multiple  wavelengths.  This  project  in  the  current  reporting  period 
involves  the  theoretical  modeling  of  photon  migration  in  tissue  using  the  cumulant  solution  to  radiative 
transfer  and  implementing  an  enhanced  3D  tomographic  image  reconstruction  algorithm.  Significant 
advances  in  both  fronts  were  made  during  the  current  reporting  period. 

Body 

The  tasks  performed  and  the  progress  made  during  th  e  current  reporting  period  include  the  theoretical 
modeling  of  photon  migration  in  tissue  using  the  cumulant  solution  to  radiative  transfer  and  implementing 
enhanced  3D  tomographic  image  reconstruction  algorithms. 

Theoretical  modeling  of  photon  migration  in  tissue 

We  extended  the  photon  transport  model  for  light  migration  in  turbid  media  based  on  a  cumulant 
approximation  to  radiative  transfer  to  a  bounded  medium  with  planar  geometries  (Task  1.1).  This 
extension  makes  the  cumulant  model  more  suitable  in  practical  applications  which  always  involve  finite 
geometries.  This  cumulant  model  of  photon  migration  was  demonstrated  to  agree  much  better  with  the 
Monte  Carlo  simulation  than  the  conventional  diffusion  approximation  at  early  times  while  both 
approximations  agree  well  with  the  Monte  Carlo  simulations  at  later  times.[l] 

We  enhanced  the  3D  tomographic  image  reconstruction  algorithms  [2,  3]  by  using  the  new  cumulant 
transport  model  (Task  1.2).  By  scanning  a  point  source  on  the  grids  of  the  input  plane  of  a  slab  and 
measuring  light  intensity  on  a  detector  array  on  the  exit  plane  of  the  slab,  a  set  of  four-dimensional  (4D) 
data  is  formed.  A  near  real-time  image  reconstructi  on  is  obtained  by  performing  a  hybrid-dual-Fourier  on 
the  4D  data  set.  [4, 5,  See  Appendix  1  and  3]  Im  provement  in  the  modeling  accuracy  was  achieved. 

To  access  the  efficacy  of  a  linear  inversion  scheme  in  image  reconstruction  of  human  breasts,  we  also 
studied  the  nonlinear  effect  of  the  multiple  passages  of  an  absorption  inhomogeneity  of  finite  size  deep 
inside  a  turbid  medium  on  optical  imaging  using  the  cumulant  solution  to  radiative  transfer  (Task  1).  We 
derived  the  analytical  nonlinear  correction  factor  wh  ich  agrees  excellently  with  the  predictions  from  the 
Monte  Carlo  simulations.  We  concluded  that  the  effect  of  the  nonlinear  multiple  passages  of  an 
absorption  site  on  optical  imaging  only  becomes  appreciable  when  the  size  of  the  inhomogeneity  reaches 
10/,  (ten  times  of  mean  free  transport  paths)  or  larger  for  human  tissues.  [6,  See  Appendix  4] 


Implementing  of  the  enhanced  3D  tomographic  image  reconstruction  algorithm 

The  non-iterative  3D  tomographic  reconstruction  algorithm  (hybrid-dual-Fourier  tomography)  was 
implemented  (Task  1.2).  The  reconstruction  uses  fast  Fourier  transforms  (FFT)  and  is  fast.  A  complete 
reconstruction  of  32  x  32  x  20  voxels  is  within  few  minutes  using  1GHz  CPU. 

We  studied  the  relation  between  the  appropriate  re  gulation  and  the  noise  presented  in  measurement  (Task 
1.4).  We  found  the  regularization  parameter  is  linked  directly  to  the  prior  signal  to  noise  ratio.  We 
provided  evidence  that  proper  modeling  of  the  noise  and  appropriate  regularization  improves  the  quality 
of  image  reconstruction.[7,  See  Appendix  5] 

We  are  working  on  including  of  multiple  wavelengths  of  probing  light  in  our  tomographic  imaging 
algorithms. 
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Key  Research  Accomplishments 

•  Extended  cumulant  solution  of  radiative  transfer  to  planar  geometries  making  it  more  suitable  for 
practical  applications  which  involve  finite  boundaries. 

•  Developed  and  enhanced  the  3D  tomographic  image  reconstruction  algorithm  by  using  the 
cumulant  transport  model. 

•  Developed  the  criterion  of  the  optimal  regulati  on  parameter  for  inverse  image  reconstruction 
related  to  the  noise  presented  in  measurements. 

•  Derived  the  nonlinear  correction  factor  of  multiple  passages  of  an  absorption  inhomogeneity  by  a 
photon  for  optical  imaging  and  provided  a  measure  of  the  efficacy  of  linear  inversion  schemes. 


Reportable  Outcomes 

Journal  Papers: 

1 .  Cai,  W.,  M.  Xu,  and  R.R.  Alfano.  Three  dimensional  radiative  transfer  tomography  for  turbid 
media,  in  IEEE  JSTQE .  2003  (to  appear  in  press) 

2.  Xu,  M.,  M.  Lax,  and  R.R.  Alfano,  Light  anomalous  diffraction  using  geometrical  path  statistics 
of  rays  and  gaussian  ray  approximation.  Opt.  Lett,  2003.  28:  p.  179-181. 

Presentations  and  Proceeding  Papers: 

3.  Xu,  M.,  W .  Cai,  and  R.R.  Alfano,  Three  dimensional  Hybrid-Dual-F ourier  tomography  in  turbid 
media  using  multiple  sources  and  multiple  detectors,  in  Third  inter-institute  workshops  on  diagnostic 
optical  imaging  and  spectroscopy:  the  clinical  adventure .  2002:  National  Institute  of  Health,  Bethesda, 
MD 

4.  Xu,  M.,  W .  Cai,  and  R.R.  Alfano.  Nonlinear  multiple  passage  effects  on  optical  imaging  of  an 
absorption  inhomogeneity  in  turbid  media,  in  European  Conference  on  Biomedical  Optics:  Photon 
migration  and  Diffuse-light  imaging.  2003 

Manuscripts: 

5.  Xu,  M.,  et  al..  Prior  information  and  noise  in  three-dimensional  optical  image  reconstruction . 

Conclusions 

The  work  carried  out  during  the  current  reporting  period  builds  on  and  affirms  some  of  our  earlier 
inferences  and  leads  to  the  following  conclusions.  First ,  the  cumulant  transport  model  provides  a  more 
accurate  model  than  the  conventional  diffusion  model  for  the  description  of  light  propagation  in  turbid 
media  such  as  human  breasts.  Second ,  the  optimal  regularization  of  image  reconstruction  depends  on  the 
noise  presented  in  the  measurements;  proper  modeling  of  the  noise  and  appropriate  regularization 
improves  the  quality  of  image  reconstruction.  Third ,  the  nonlinear  effect  of  the  multiple  passages  of  an 
absorption  site  by  a  photon  on  optical  imaging  only  becomes  appreciable  when  the  size  of  the 
inhomogeneity  reaches  10/,  or  larger  in  human  tissues.  Fourth,,  the  theoretical  formalism  and  computer 

algorithm  for  3D  tomographic  image  reconstruction  s  hows  (with  simulated  data)  the  potential  to  provide 
fast  3D  images  of  the  scattering  and  absorption  objects  at  various  depths  in  turbid  media. 
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Abstract 

The  photon  distribution,  as  a  function  of  position,  angle  and  time,  is  computed  using  the 
analytical  cumulant  solution  of  the  Boltzmann  radiative  transfer  equation  (RTE).  A  linear  forward 
model  for  light  propagation  in  turbid  media  for  three  dimensional  (3D)  optical  tomography  is 
formed  based  on  this  solution.  The  model  can  be  used  with  time  resolved,  CW,  and  frequency- 
domain  measurements  in  parallel  geometries.  This  cumulant  forward  model  (CFM)  is  more 
accurate  than  that  based  on  the  diffusion  approximation  of  RTE.  An  inverse  algorithm  that 
incorporates  this  CFM  is  developed,  based  on  a  fast  3D  hybrid-dual -Fourier  tomographic  approach 
using  multiple  detectors  and  multiple  sources  in  parallel  geometries.  The  inverse  algorithm  can 
produce  a  3D  image  of  a  turbid  medium  with  more  than  20,000  voxels  in  1-2  minutes  using  a 
personal  computer.  A  3D  image  reconstructed  from  simulated  data  is  presented. 

Subject  terms:  photon  migration;  radiative  transfer  equation;  forward  model; 

absorption  and  scattering;  optical  tomography;  inverse  algorithm. 


Accepted  by  IEEE  special  issue  JSTQE  (Laser  in  medicine  and  biology)  #1277 

1.  Introduction 

Over  the  past  decade,  optical  tomography  has  been  investigated  as  a  nonin vasive  imaging  method  that  uses 
non-ionizing  near-infrared  (NIR)  light  to  obtain  images  of  the  interior  of  the  breast.  Unlike  X-ray,  which  is 
attenuated  through  media  by  ionizing  the  electrons  at  inner-orbits  of  atoms,  NIR  light  uses  the  vibrational  overtones 
for  different  molecular  components  in  the  structures  of  tumor.  NIR  light  may  be  used  to  create  image  based  on  the 
molecular  change,  which  may  be  used  to  improve  sensitivity  and  specificity  in  the  early  diagnostics  of  breast  cancer. 
Breast  tissues  scatter  light  strongly,  and  blur  the  direct  shadow  image  of  a  tumor.  A  technique,  known  as  inverse 
image  reconstruction,  has  been  investigated  to  overcome  the  problem  of  multiple  scattering.  Some  obstacles  in  the 
development  of  optical  tomography  are  inaccuracy  of  the  commonly  used  diffusion  forward  model,  and  lack  of  a 
fast  inverse  algorithm  able  to  realize  a  three  dimensional  (3D)  image  reconstruction  of  a  breast  for  clinical  use. 

One  critical  issue  is  the  forward  model,  which  should  correctly  simulate  photon  propagation  in  the  medium. 
The  most  commonly  used  forward  models  were  built  based  on  solution  of  the  diffusion  equation,  which  is  the  lowest 
approximation  of  the  radiative  transfer  equation  (RTE).[l-5]  The  forward  models  based  on  the  diffusion 
approximation  (DA)  give  a  large  error  when  the  distance,  d ,  between  a  voxel  and  a  source  is  small.  Furthermore,  the 
photon  distribution  still  maintains  a  strong  anisotropy  in  a  deeper  region  away  from  a  source,  which  will  be  shown 
later  in  this  paper.  Unfortunately,  contributions  from  near  surface  voxels  to  measured  signals  are  often  larger  than 
contributions  from  the  voxels  deep  inside  the  medium.  Inaccuracy  of  DA  based  forward  model  may  lead  to  a  failure 
in  image  reconstruction,  especially  for  small  hidden  objects  deep  inside  the  medium.  The  total  weight  matrix  should 
be  inverted.  The  large  elements  in  the  matrix,  which  play  a  more  important  role  in  inversion,  are  evaluated 
incorrectly  in  DA  models.  The  shortcoming  of  DA  is  well  recognized,  but  it  is  still  broadly  applied  due  to  the 
difficulty  in  directly  solving  the  radiative  transfer  equation.  Hielscher  et  al  [6]  and  Vihunen  et  al  [7]  developed 
numerical  solutions  of  RTE  for  optical  tomography. 
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Recently,  we  have  developed  an  analytical  solution  of  RTE,  based  on  cumulant  expansion,  in  an  infinite 
uniform  medium  with  an  arbitrary  phase  function.  [8,9]  It  provides  an  explicit  analytical  expression  for  photon 
distribution  function  /( r,  s,  t\  as  a  function  of  position  r,  direction  of  light  s,  and  time  t  The  mean  position  and  the 
half-width  at  half  maximum  height  (HWHM)  of  the  distribution  are  always  exact.  In  this  paper,  the  linear  forward 
model  based  on  the  cumulant  solution  is  described.  This  CFM  may  used  with  time-resolved,  CW,  and  frequency- 
domain  data,  which  are  much  more  accurate  than  the  DA  models. 

To  obtain  a  3D  image  one  needs  to  investigate  the  inverse  algorithms.  For  clinical  applications,  this  requires  an 
inversion  technique,  that  is  computationally  fast,  and  stable  in  the  presence  of  measurement  noise.  Recent 
algorithms  to  solve  the  inverse  problem  include  Newton’s  least-square-based  methods  and  gradient-descent 
methods. [1-5]  These  approaches  use  an  iterative  procedure,  which  requires  a  long  computation  time  to  solve  a  3D 
inverse  problem  with  large  unknowns  (the  number  of  unknowns  is  the  number  of  voxels).  Furthermore,  the  iterative 
methods  can  not  ensure  that  the  result  arrives  at  a  “global  minimum”,  and  does  not  converge  to  a  “local  minimum”, 
which  is  not  a  true  image  of  the  medium.  The  application  of  Fourier  transform,  which  has  been  called  ’’diffraction 
tomography”,  can  greatly  reduce  computation  time.  Matson  et  al  [10]  and  Li  et  al  [11]  have  developed  the 
diffraction  optical  tomographic  methods  to  realize  fast  image  reconstruction.  However,  their  algorithms  are  limited 
to  the  use  of  a  single  light  source  with  a  2D  plane  of  detectors.  This  type  of  experimental  setup  acquires  only  a  set  of 
2D  data  using  CW  or  frequency-modulated  light,  that  is  not  enough  for  a  3D  image  reconstruction.  Recently, 
Schotland  and  Markel  developed  inverse  inversion  algorithms  using  diffusion  tomography  [12,13,14]  based  on  the 
analytical  form  of  the  Green’s  function  of  frequency-domain  diffusive  waves,  and  point-like  absorbers  and 
scatterers.  Data  obtained  by  multiple  sources  with  multiple  detectors  in  parallel  slab  geometiy  are  used  in  these 
approaches. 

A  fast  hybrid-dual-Fourier  (HDF)  algorithm,  which  uses  multiple  sources  and  multiple  detectors  in  parallel  slab 
geometry,  is  described  in  this  paper  for  reconstruction  of  a  3D  image  of  an  inhomogeneous  medium.  This  approach 
uses  a  general  2D  translation  invariance  of  the  Green’s  function  in  a  homogeneous  background  slab  medium, 
suitable  for  forward  models  based  on  solution  of  RTE,  and  various  other  forward  models,  in  CW,  frequency-domain, 
and  time-resolved  measurements.  This  inverse  algorithm  runs  fast.  It  is  shown  that  a  3D  image  of  a  turbid  medium 
(for  example,  divided  into  32x32x20=20480  voxels)  can  be  reconstructed  in  1-2  minutes  using  a  personal  computer. 
This  algorithm  can  produce  stable  images  in  presence  of  relatively  strong  noises. 

The  forward  model  and  the  inverse  algorithm  discussed  below  can  also  be  applied  for  image  reconstruction  in  a 
cloudy  environment  for  military  use. 

This  paper  is  organized  as  follows:  Section  2  presents  the  analytical  solution  of  RTE,  based  on  a  cumulant 
expansion,  in  an  infinite  uniform  medium  and  shows  the  photon  distribution  function  computed  using  the  cumulant 
analytical  solution.  Section  3  describes  the  forward  models  based  on  the  analytical  solution  of  RTE,  considering  the 
slab  geometry,  and  a  weak  heterogeneity  using  a  perturbative  method.  Section  4  describes  the  hybrid-dual-Fourier 
inverse  algorithm  for  a  reconstruction  of  a  3D  image  of  an  inhomogeneous  medium.  The  3D  image  using  this 
algorithm  is  shown.  A  discussion  is  presented  in  Section  5. 

2.  Analytical  cumulant  solution  of  RTE 

The  photon  propagation  in  a  medium  is  described  by  the  photon  distribution  function,  /( r,  s,  t\  as  a  function  of 
time  t,  position  r,  and  direction  s.  The  mathematical  equation  governing  photon  propagation  is  the  well-known 
radiative  transfer  equation: 

dl  (r,  s,  0  /  dt  +  cs  •  V  r /(r,  s,  t )  +  fia  (r)/(r,  s,0  =  /^(r  )|  P(  s,  s’ ,  r)[/(r,  s’ ,  /)  -  7(r,  s,  /)]<*’ 

+  S(r  -r0)£(s-s0 )5(t  -  0)  (1) 

where  the  fundamental  parameters  are  the  scattering  rate  ps(r)  =  cpcrs,  the  absorption  rate  p„(r)  =  cpoa,  and  the 
differential  angular  scattering  rate  ps(r)P(s,  s',  r),  where  o„  and  os  are  the  absorption  and  scattering  cross  sections 
respectively,  p  is  density  of  scatterers,  and  c  is  the  speed  of  light  in  the  medium.  In  a  uniform  infinite  medium,  these 
parameters  are  position  independent. 

When  the  phase  function  depends  only  on  the  scattering  angle,  we  can  expand  the  phase  function  in  Legendre 
polynomials  with  constant  coefficients, 

P( S,s’)  = - cos(s-s’)]  (2) 

4  7t 

Recently,  we  have  developed  a  new  approach  to  obtain  an  analytical  solution  of  RTE,  based  on  a  cumulant 
expansion,  in  an  infinite  uniform  medium,  with  an  arbitrary  phase  function  P( s,  s').  [8,9] 
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We  briefly  review  the  concept  of  “cumulant”  in  a  ID  case.  Consider  a  random  variable  x,  with  a  probability 
distribution  function  fix).  Instead  of  using  fix)  to  describe  the  distribution,  we  define  the  «th  moment  of  x , 

<xn>=  ^xnf{x)dx>  and  correspondingly  the  wth  cumulant  <x”>c  defined  by 


OO  00 

exp(y <  x  >c  (it)  / n\)  =<  exp(zYx)  >-^<xn  >(it)n  / n\.  The  first  cumulant  <x>c  is  the  mean 

n- 1  *= 0 

position  of  x.  The  second  cumulant  <x2>c  represents  the  HWHM  of  the  distribution.  The  higher  cumulants  are  related 
to  the  detailed  shape  of  the  distribution.  For  example,  <x3>c  describes  the  skewness  or  asymmetry  of  the  distribution, 
and  <x4>c  describes  the  “kurtosis”  of  the  distribution,  that  is  the  extent  to  which  it  differs  from  the  standard  bell 
shape  associated  with  the  normal  distribution  function.  The  cumulants,  hence,  describe  the  distribution  in  an 
intrinsic  way  by  subtracting  off  the  effects  of  all  lower  order  moments.  In  3D  case,  the  first  cumulant  has  3 
components,  the  second  cumulant  has  6  components,  and  so  on. 


We  derived  an  explicit  algebraic  expression  of  spatial  cumulants  at  any  angle  and  any  time  that  is  exact  up  to  an 
arbitrarily  high  order  n .  [9]  This  means  the  distribution  function  /(r,  s,  t)  can  be  computed  to  any  desired  accuracy. 
At  the  second  order,  n  =  2,  an  analytic,  hence,  useful  explicit  expression  for  distribution  function  7( r,  s,  /)  is 
obtained.  [8]  This  distribution  is  Gaussian  in  position,  which  is  accurate  at  later  times,  but  only  provides  the  exact 
mean  position  and  the  exact  HWHM  at  early  times.  A  weakness  of  the  second  order  cumulant  solution  is  that 
photons  at  the  front  edge  of  Gaussian  distribution  travel  faster  than  light  speed,  thus  violate  causality,  though  to  a 
much  less  extent  than  that  in  the  diffusion  approximation. 

Fig.  1  compares  7(r,  s,  t)  obtained  from  the  analytical  cumulant  solution  and  the  Monte  Cairo  simulation.  In 
order  to  reduce  the  statistical  deviation  to  an  acceptable  level,  109  events  are  counted  in  the  Monte  Carlo  simulation. 
The  figure  shows  that  die  solid  curve  (the  10th  order  cumulant  solution)  is  located  in  the  middle  of  data  obtained  by 
the  Monte  Carlo  simulation.  The  solution  for  CW  case  can  be  obtained  by  an  integration  of /(r,  s,  t)  over  time  t.  It  is 
shown  that  even  second  order  cumulant  solution  (the  dotted  curve)  can  provide  an  accurate  CW  solution,  because 
this  solution  ensures  that  the  mean  position  and  the  HWHM  of  distribution  are  always  exact. 

The  plots  in  Fig.  1  indicates  that  a  strong  anisotropic  angular  distribution  still  exists  at  z  -  6  ltr  (/^  is  the 
transport  mean  free  path)  from  the  source.  The  diffusion  approximation  is  only  valid  when  the  angular  distribution  is 
nearly  isotropic.  The  dominate  5  wave  distribution  N( r,  t)/4n  computed  using  the  diffusion  model  (the  thick  dotted 
curve)  has  a  large  discrepancy  with  the  Monte  Carlo  result. 

The  second  order  analytical  cumulant  solution  is  given  by  [8] 


7(r,s,0  = 


F(s,s0,  t) 


1 


(4 7Z)V2  (det B) 


Ml 


1 

exp  — 

L  4 


)a(r-rc). 


(3) 


where 


2/  +  1 

F(s,s0,O  =  exp  (-//„/)£,  — — exp(-g/f)/>,[cos(s  •  s0 )] 

4  7U 

2/  +  1  *  (4) 

=  exp(-/V)E/— -  expC-S/OE  Ylm  (s  )Ylm  (s0 ) 

4  7t  m 

In  Eq.  (4),  g}  =  ps[l  -  aj( 2/+1)],  F/m(05(|))=(— l)m[(/— m)I/((/+A?r)!]I/2/>/('”)(cos0)exp(//?2(|>),  where  P/m)( cos0)  is  the 
associated  Legendre  function,  and  Y,m  (s)  are  spherical  harmonics  normalized  to  47t/  (2/+1). 

In  Eq.  (3),  the  mean  position  of  the  distribution  (first  cumulant),  when  the  source  is  located  at  r0  =  0  and  the 
incident  direction  is  along  z,  is  given  by: 

<  (s>  0  =  GZ,  A,p,  (cos  0)[(l  + 1)/ (g,  -  gM  )  +  If  (g,  -  g,_, )] 

r;(s,t)  =  G'£lAlP")(cos0)cos0[f(g, -gM)-/(g,  -#/+l)] 

where  G  =  c  exp  (-|V)/F(s,  s0,  t ),  A,  =  (l/47t)exp(-g,r),  and 

/(g)  =  [exp(gO-l]/g 

r,  "  is  obtained  by  replacing  cos<(>  in  Eq.  (5.2)  by  sin<j>- 
The  HWHM  (second  cumulant)  is  expressed  as 


(5.1) 

(5.2) 

(6) 
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Bap  (s»  0  =  cGAag  -  r'rcg  /  2  , 


Aa  =  XiAfii cos0)  EjV  +-£-£«  +^U<4>  .  CM) 

A„.„  =Z,|  «(cose|-^£«'>- +  t2£f>  +  y^v 

±z^^a>(C0S#)CM(^^£'<'>+^£",-^£”>-^^,>]  •  <8'2> 


where  (+)  corresponds  to  A**  and  (-)  corresponds  to  A>: 

A.„.  =  1/  4  4^ (2)  (cos  <9)  sm(2^)|" — E\X)  +  - 


>+— £,(2) - —  £/3) - —  £<4>1,  (8.3) 

\  1  1  '*W  1  ‘  >vy 


A,z  =I/|4^(,)(cos^)cos(^) 


— £/4) 

/  +  3  J 


2 . L2/-1  '  2/  +  3  '  2/  —  1  2/  +  3  J  ’  v  ' 

A.  =  I,|^f,)(cos^)cos(^)  +^I£,(3)  +^£/4)j  (8.4) 

is  obtained  by  replacing  cosc|>  in  Eq.  (8.4)  by  sin<{).  In  Eqs.  (8.1-  8.4)  £/(l_4)  are  given  by: 

E?]  =  U(g,-g,-2)-f{g,  -g,_i)]/(gM -g,-z)  ,  (9.1) 

EI2)  =  [f{g,  - g/+2 ) -  f{g i  - gM )]  fgM  - g I +2 )  ,  (9.2) 

E?]  =  Uig,  ~  g,-, )  -  f  !{gi  ~  g,_, )  ,  (9.3) 

E)*  =lf{g-gM)-t]l{g,-gM)  .  (9.4) 

Figs.  2(a)  and  2(b)  show  the  light  distribution  as  a  function  of  time  at  different  receiving  angles  in  an  infinite 
uniform  medium,  computed  by  the  second  cumulant  solution,  where  detector  is  located,  separately,  at  5  ltr  (Fig.  2a) 
and  15  ltr  (Fig.  2b)  from  the  source  in  the  incident  direction  of  the  source.  Fig.  2  shows  the  existence  of  the  strong 
anisotropy  of  the  light  distribution  at  5  ltr  from  the  source  and  the  modest  anisotropy  at  a  distance  of  15  /*  .  These 
types  of  distributions  have  been  demonstrated  by  time-resolved  experiments.  [15] 

One  advantage  of  using  the  above  analytical  solution  of  RTE  is  that  the  distribution  function  can  be  computed 
very  fast.  The  associated  Legendre  functions  can  be  accurately  computed  using  recurrence  relations.  It  takes  only  a 
minute  to  compute  105  data  of  7(r,  s,  0  on  a  personal  computer. 

The  corresponding  solution  in  the  frequency-domain  /( r,  s,  co)  can  be  obtained  by  making  a  Fourier  transform 
| dt  exp(— iO)t)I(r,S,t)  .  The  CW  solution  is  obtained  by  taking  0. 

The  photon  density  N( r,  /)  of  the  second  cumulant  solution  is  given  by 


^  1  1  (z-Rrf 

N(ryt)  = - - -pr- - exp - exp  - 

(4 iiDzzctY2  InDxxCt  [  4Dzzct  J 

with  the  mean  position  Rz  —  c[l  —  exp(— g-^)]/^  . 

The  corresponding  time-dependent  diffusion  coefficients  are: 


(x2  +y2) 
4  Dvxt 


exp  (~Ma‘)  * 


f|-  +  3g,  g2  [i_exp(_gi0]  + 


3Mgl  g|(gl-g2) 


g2(gl-g2> 


[1  -exp(-g20] - j[\ -expf-g,?)]2  [  (12) 


D^=Dyy=YA~  +~TT1 - :  [  1  -  exP(-gi; )]  +  — T— - -[1  — exp(— g-,r))l  ■  (13) 

3?Ui  gr(gi-g2)  g2(gi-g2>  ■  j 

As  shown  in  Eqs.  (1 1)  -  (13).  the  mean  position  of  the  distribution  is  moving,  and  the  diffusion  coefficients  are 
time  dependent.  At  t  — *  0,  the  mean  position  of  the  photon  density  moves  along  z  direction  with  speed  c,  and  the 
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diffusion  coefficients  tend  to  zero,  this  result  presents  a  clear  picture  of  near  ballistic  motion.  As  time  increases,  the 
mean  position  motion  slows  down,  and  the  diffusion  coefficients  increase  from  zero.  This  stage  of  photon  migration 
is  often  called  a  snakelike  mode.  At  long  time,  Eq.  (10)  tends  to  the  center-moved  (1  ltr)  diffusion  model  with  the 
diffusion  coefficient  IJl. 

3.  Forward  model  based  on  the  cumulant  solution  of  RTE 

The  linear  forward  models  for  scattering  media  are  built  in  following  three  steps:  (1)  computation  of  a 
background  Green’s  function  in  an  infinite  uniform  medium;  (2)  extension  of  this  Green’s  function  to  slab 
geometry;  and  (3)  computation  of  the  weight  function  using  a  perturbative  method.  These  steps  have  been  applied  in 
building  the  linear  forward  models  under  DA.  [2]  We  use  these  steps  as  well,  but  our  approach  is  based  on  the 
cumulant  solution  of  RTE,  rather  than  the  solution  of  the  diffusion  equation. 

We  use  the  second  order  cumulant  solution  for  computing  a  background  Green’s  function  in  an  infinite 
uniform  medium,  since  it  is  easy  to  use  the  explicit  expressions  in  Eq.  (3)  -  Eq.  (9),  that  avoid  complicated 
computations  of  higher  order  cumulants.  The  second  order  cumulant  solution  is  accurate  at  later  times,  but  only 
provides  the  correct  mean  position  and  the  correct  HWHM  at  early  times.  We  notice  that  the  width  of  the 
distribution  at  early  times  could  be  smaller  than  the  size  of  a  voxel,  the  average  over  the  distributions  at  different 
points  in  a  voxel  smears  the  detail  shape  of  the  distribution.  In  the  CW  or  frequency-domain  cases,  the  shape  of  the 
distribution  is  further  smeared  by  integration  over  time  t  Therefore  the  second  order  cumulant  solution  can  be  a 
reasonable  approximation  in  building  forward  models  based  on  the  RTE. 

Since  a  detector  usually  collects  emergent  light  within  a  wide  range  of  angle  of  different  directions,  it  is 
convenient  to  compute  the  Green’s  function  related  to  a  detector  using  photon  density  N(rd,  t)  [Eqs.  (10)-(13)], 
where  ra  is  the  position  of  detector. 

It  is  essential  to  include  the  boundary  effect  in  the  solution  of  the  RTE  when  photons  are  injected  into  and 
spread  out  from  a  finite  sized  medium.  A  proper  extension  of  the  cumulant  solution  to  slab  geometry  is  an  essential 
step  for  building  a  forward  model. 

A  boundary  condition  is  applied  based  on  the  following  physical  consideration.  At  early  times,  the  center  of 
photon  distribution  injected  into  medium,  moves  forward  into  medium.  Then  the  distribution  spreads  out  from  the 
moving  center  with  diffusion  coefficients  that  gradually  increase  from  zero.  At  early  times,  the  number  of  photons 
leaking  out  of  the  boundary  is  negligible  compared  to  the  total  number  of  the  incident  photons.  The  boundary 
condition  plays  a  role  at  later  times,  when  there  are  many  photons  leaking  out  of  the  boundary. 

The  approach  known  as  an  approximate  ’’extrapolated"  boundary  condition  [16],  extrapolates  the  boundary  by 
a  distance  £  =  a  ltn  the  extrapolation  length,  beyond  the  real  boundaries  with  a  ~  0.7,  at  which  the  photon  density 
vanishes. 

To  apply  this  boundary  condition  for  the  cumulant  solution  in  a  semi-infinite  geometry,  a.  virtual  negative 
source,  Sv,  is  added  to  the  original  source,  S,  as  shown  in  Fig.  3.  During  the  early  period,  the  solution  of  the  RTE  in 
an  infinite  uniform  medium  automatically  satisfied  the  boundary  condition  because  the  density  is  near  zero  at  the 
boundary,  and  the  virtual  source  does  not  play  a  role.  After  a  time  of  approximately  4  IJc ,  the  center  of  photon 
density,  C,  has  moved  and  stopped  at  a  position  1  from  the  original  source  S  and  the  center  from  virtual  source, 
Cv,  has  moved  in  a  similar  way.  Then,  the  arrangement  shown  in  Fig.  3,  produces  a  cancellation  of  contributions  to 
the  photon  density  from  the  original  source  and  the  virtual  source  on  the  extrapolated  boundary. 

Fig.  4  shows  that  the  time-resolved  backscattered  photon  distribution  in  a  semi-infinite  medium  on  the  z  =  0 
surface,  with  the  source-detector  distance  1  ltn  obtained  using  the  second-order  cumulant  approximation  and  the 
extrapolated  boundary  condition,  which  agrees  with  the  Monte-Carlo  simulation  much  better  than  that  of  the  DA. 

For  extending  to  the  slab  geometry,  adding  a  series  of  pairs  of  virtual  "image”  sources  at  both  sides  of  slab  is  a 
good  approximation  for  satisfaction  of  the  extrapolated  boundary  conditions  on  both  sides  of  a  slab.  [17] 

The  heterogeneous  structure  of  a  highly  scattering  turbid  medium  can  be  characterized  by  the  following 
optical  parameters:  the  scattering  rate  ps(r),  the  absorption  rate  pfl(r),  and  the  differential  angular  scattering  rate 
ps(r)P(s,  s’,  r). 

A  perturbation  method  is  used  which  takes  the  photon  distribution  function  in  a  uniform  background  slab 
medium  as  the  zero-order  approximation.  The  change  of  the  photon  distribution  function  originates  from  the  change 
of  optical  parameters  compared  to  that  in  the  uniform  background  slab  medium.  The  change  of  scattering  and 
absorption  parameters  are  defined  as  follows: 

Au.(r)  =  |is(r)  -  n/0) , 

A|i„(r)  =  n?(r)  -  |io(0) ,  (14) 

A[|.isP](s,  s',  r)  =  |is(r)P(s,  s’,  r),  -  ii/V0)(s,  s'). 
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where  the  quantities  with  super  index  (0)  are  the  optical  parameters  in  a  uniform  background  slab  medium.  By 
expanding  A[psP](s,  s ,  r)  in  Legendre  polynomials,  we  obtain 

A[yUsP](s,s',r)  =  (r>/0)  +/u/0)Aa/(r)]P;[cos(s -s')]  ,  (15) 

with  Aa0(r)  =  0,  since  a0  always  equals  to  1.  The  physical  meaning  is  that  the  scattering  parameters  have  no  effect  on 
the  s  (/  =  0)  component. 

Making  a  perturbation  expansion  of  Eq.  (1)  to  the  first-order  Bom  approximation,  the  change  in  the  photon 
distribution  is  given  by 

&I(rd,sd,t\Ts,ss)=  jdt'  j<*  Jds7(0)  (rd  ,sd  ,t  - 1'\  r,  s') 

{ jA[psP](s,s’,r)/(0)(r,s,P  |  rs,ss)ds - (r)  +  A//a (r)]/(0) (r,s V'|  rs ,ss )}  ’  ' 

where  A I  (rAsA/|rs,ss)  is  the  change  in  the  light  intensity  received  by  a  detector  located  at  rd,  along  the  direction  sdy 
and  at  time  t,  which  is  injected  from  a  source  located  at  rs,  along  a  direction  of  ss,  at  time  t  =  0.  "Change”  refers  to 
the  difference  in  intensity  compared  to  that  received  by  the  same  detector,  from  the  same  source,  when  light  passes 
through  a  uniform  background  slab  medium.  The  term/0)  (r2,s2,/|ri,S|)  is  the  intensity  of  light,  calculated  using  the 
cumulant  solution  of  RTE,  at  r2  along  the  direction  s2  and  at  time  t,  when  light  is  injected  from  a  position  i*i  along  a 
direction  of  S|  at  time  t  ~  0  migrating  in  a  uniform  background  slab  medium. 

The  background  Green’s  functions  in  Eq.  (16),  obtained  by  cumulant  solution,  are  expanded  in  spherical 
harmonics: 

/<0)(r,M’l  r„s,)  =  2^(r»r,,s„Orfa(*). 

l,m 

I(0\rd,sdd-t'\r,s)  =  YCL(r>rd’sd>t-tlYL(s)  ■  (17) 

I,m 

The  spherical  transform  is  performed  using  a  fast  Fourier  transform  for  the  integral  over  <|>,  and  a  Clenshaw-Curtis 
quadrature  for  the  integral  over  G. 

Using  the  orthogonality  relation  of  the  spherical  function  and  the  addition  theorem: 
(s0  —  [cos(s  •  S r)] ,  the  analytical  integration  over  s  and  sr  in  Eq.  (16)  can  be  performed.  For 

m 

time  resolved  data,  the  contribution  from  an  absorbing  object  located  at  rk  is  given  by 


U(rd,sd,rs,ssJ\rk)  =  -Ava(rk)SVkldt'£J— 

1=0  +  f )  m 

where  SVk  is  the  volume  of  kih  voxel,  and  L  is  the  cut-off  value  in  the  Legendre  expansion  in  Eq.  (18).  The 
contribution  from  a  scattering  object  located  at  rk  is  given  by 

AI(rd,sd,rs,ss,t\rk)  = 


■dl\  I A//  (rk )  1-^i— 
*  tf(2l  +  Y)  1  Al  21  +  1 


For  Frequency  domain  (or  CW)  data,  the  contribution  from  an  absorbing  object  located  at  r*  is  given  by 

A/(rrf,s(/,ri,ss,o|ri)  =  -A na (rt  )SVk £  — —  ^Alm (r, ,r.,ss, o)cbn (r4 ,  r , ,  s d,o>),  (20) 

/=0  1  m 

and  the  contribution  from  a  scattering  object  located  at  r*  is  given  by 

M(rd,sd,rs,s5,a)\rk)  = 

*  /  Nf,  a<i0)  rot  A«/(ri  )lv-«  ,  .  (2i; 

_^Z^7TT  AVs(rk)  l-TT—  -mV  —  ■  .  YjAiVrk^s^s,0))cln,{rk,rd,sd,(o). 
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Comparing  Eqs.  (18) — (21)  with  the  corresponding  weight  function  commonly  used  in  the  diffusion 
approximation,  [1,2]  only  s  wave  (/=0)  for  absorptive  objects,  and  only  p  wave  (M)  for  scattering  objects  are 
considered  in  the  diffusion  forward  models.  Besides,  even  for  s  wave  and  p  wave,  the  diffusive  solution  is  incorrect 
when  voxels  are  located  near  the  source,  as  discussed  before. 

Above  formulae  allow  simulating  the  background  Green’s  function  and  the  change  of  optical  parameters  in 
detail.  They  are  also  applicable  to  the  cases  where  only  a  few  parameters  of  the  medium  are  known,  similar  to  that 
for  the  diffusion  forward  model.  When  only  and  g-factor  for  an  uniform  background  medium  are  given, 

the  Henyey-Greenstein  phase  function  [18]  is  widely  adopted  as  an  approximate  phase  function: 

nmse)  ■  =i^(2l+ i)g'p>(cosff)  ■  (22> 

Although  Eq.  (22)  uses  a  single  parameter,  g-factor,  to  describe  a  phase  function,  this  description  is  much  better  than 
that  used  in  the  DA,  which  implies  a  phase  function  linear  in  cos0  . 

If  A a,  (r)  in  Eq.  (21),  which  represent  the  change  of  the  phase  function,  is  not  considered,  two  optical 
parameters  being  imaged  are  Ap„(r)  and  Aps(r).  The  reduced  scattering  coefficient  Aps(l-ait0)/3)  is  directly  related 
to  AD  (change  of  the  diffusion  coefficient)  used  in  the  DA  models.  The  CFM,  hence,  can  be  applied  to  the 
experimental  data  in  a  similar  fashion  as  that  for  the  DA  models,  to  obtain  images  of  the  optical  parameters.  In  the 
CFM,  however,  all  contributions  from  higher  spherical  waves  are  properly  included. 

■Hie  most  time  consuming  part  in  computation  of  CFM  using  the  above  formulae  is  to  build  a  database  of  A,m 
and  C* tm.  Once  it  is  built  for  a  uniform  background  medium,  the  database  can  be  applied  for  imaging  of  various 
heterogeneity  cases.  In  parallel  geometry,  A,„  is  a  function  of  (xk-xs>  yk-ys)  due  to  the  2D  translation  invariance. 
Since  position  of  source  zs  and  incident  direction  ss  are  fixed,  only  a  3D  (xt-xs,  yk-ys,  zk)  database  is  required.  When 
ss  is  taken  along  z  direction  (light  is  injected  perpendicular  to  surface),  the  scale  of  database  is  reduced  to  2D  due  to 
the  z  axis  symmetry.  Photons  from  different  directions  in  a  wide  solid  angle  are  received  by  a  detector,  as  discussed 
before,  photon  density  N(rd-  r,  s,  t)  is  used  for  computing  the  Green’s  function  associated  with  detectors,  which  is 
independent  of  Sj ,  and  C’,„,  can  be  computed  much  easily.  The  database  can  be  built  in  a  reasonable  computation 
time  because  the  distribution  function  /0)(r2,s2,<|ri,S|)  can  be  rapidly  calculated  using  the  analytical  expressions. 

4.  Fast  3D  hybrid  dual  Fourier  (HDF)  inverse  algorithm 

We  now  outline  an  inverse  algorithm  to  quickly  reconstruct  image  of  a  medium  from  acquired  measurements 
using  the  above  CFM.  The  above  model,  neglecting  the  irrelevant  parameters,  can  be  briefly  written  as 

Y(ri,rs,zd,zs)  =  fcffdzW(rd-r,rs-r,z,zd,zs)X(r,z),  (23) 


where  R  =  (7,  z)  is  the  position  of  a  voxel  inside  turbid  medium;  f  is  (x,  y)  coordinates;  Rs  =  (?s,zs)is  the 
position  of  a  source;  and  R  d  =  (?d ,  zd  )  is  the  position  of  a  detector.  In  Eq.  (23),  Y  (?d  ,\,zd,zs)  is  the  measured 
change  in  light  intensity  received  by  a  detector  at  Rd  from  a  point  source  at  Rs .  X(r ,z)  is  the  change  of  the 


optical  parameters  inside  turbid  medium.  The  weight  function  W (rd  —  7,  ?s  —  7 ,  Z,  zd ,  zs )  is  a  function  of  7d  -  7 

and  7S  —  7  on  (x,  y)  plane,  because  of  parallel  geometry,  assuming  an  infinite  sized  area,  and  the  2D  translation 

invariance  of  the  Green’s  function  in  a  background  homogeneous  slab.  Here,  the  special  form  of  the  weight  function 
is  not  relevant;  the  weight  function  can  be  calculated  by  the  CFM  or  the  DA  models,  using  with  CW,  frequency,  or 
time-resolved  data.  This  approach  is  general  and  can  also  be  used  for  inverse  problems  of  non-optical  measurements 
in  parallel  geometries. 

A  light  source  scans  through  a  2D  array.  Transmitted  or  backscattered  light  signals  emerging  from  the 
medium  are  detected  using  a  2D  array  of  detectors,  such  as  a  CCD  camera  (or  time-gated  CCD  camera  in  the  time 
resolved  case).  Each  illumination  of  the  light  source  provides  a  set  of  2D  data  on  the  two-dimensional  detector  array. 
For  CW  or  frequency-modulated  light  source,  this  arrangement  can  produce  a  set  of  2Dx2D  =  4D  data  in  a  relatively 
short  acquisition  time,  because  a  CCD  camera  produces  2D  data  of  the  detectors  at  different  positions 
simultaneously.  When  time-resolved  or  modulation  at  multiple  frequencies  are  applied,  a  set  of  5D  data  can  be 
acquired.  The  inverse  problems  of  3D  imaging,  hence,  are  over-determined,  which  is  necessary  for  obtaining  an 
accurate  3D  image. 

When  the  translation  invariance  is  satisfied,  the  Fourier  transform  approach  is  a  powerful  technique  to  achieve 
a  fast  inversion.  In  the  Fourier  space,  the  convolution  of  W  and  X  becomes  a  product  of  W  and  X,  and  the  weight 

{ 
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matrix  W  becomes  diagonal.  Hence,  inversion  can  be  performed  much  faster.  Using  this  concept  in  the  case  of 
multiple  sources  and  multiple  detectors  in  parallel  geometries  a  dual  2D  Fourier  transform  is 

performed  on  Eq.  (23),  to  obtain 
iq  r  zq  ,  r , 

jdi ^dr^e  sse  d  d ,z^)  = 

iq  (r  -  r)  zq  ,  (r '  -  r)  z(q  +  q  , )? 

\dz\di\d(j^  -rW(?d  -r)e  s  s  e  d  d  IF(rg  - r,rd -r,z^,zd>z)e  s  d  X(r,z) 

which  leads  to 

Y{^iA,^d,zs)  =  J^(qd,qs,z,zrf>zJX(qd+qs,z),  (24) 

where  F,  1,  and  W  are  change  in  light  intensity,  change  in  optical  parameters,  and  the  weight  function  in  the 
Fourier  space  respectively. 

A  similar  form  of  this  dual  Fourier  transform  has  been  derived  by  Markel  and  Schotland  [1 3,14]  in  a  frequency- 
domain  diffusion  model. 

Eq.  (24)  seems  difficult  to  be  used  for  performing  the  inverse  reconstruction  because  of  the  argument 
mismatch  (qd  -f  qs)  in X  and  (qs ,  qd  )  in  Y  and  W .  This  difficulty  occurs  because  the  weight  function  in  Eq. 

(23)  is  related  to  three  positions:  rd ,  rs ,  and  ?  .  To  remove  this  complexity,  the  following  linear  hybrid  transform  is 
introduced: 


u  =  qd  +qs 

v  =  qd-qs 

This  results  in  HDF  formula: 


Y  (u,  v,zd,zs)  =  jdzW{ u,  v,  z,  zd ,  zs  )X  (u,  z) 


(25) 

(26) 


where  Y  ,  X ,  and  W  are,  respectively,  Y ,  X ,  and  W  as  functions  of  u  and  v  . 

While  Eq.  (25)  is  a  relatively  simple  expression,  it  is  essential  to  properly  realize  this  hybrid  transform  in 
discrete  lattices  of  the  Fourier  space.  A  procedure  to  quickly  perform  this  transform  from  (qd,  qs)  coordinates  to  new 
(u,  v)  coordinates,  separately,  for  x  and  y  components,  is  explained  in  Fig.  5  using  an  example  of  a  6x6  lattice.  The 
maximum  value  of  u  is  taken  as  the  maximum  value  of  qd  or  qSj  not  the  maximum  value  of  qd  +  qs.  The  periodic 

property  of  lattice  in  the  Fourier  space  is  used,  for  example,  Y  (u=2,  v=4)  =  Y  (qd=3,  q=5).  This  procedure  builds  a 
one-to-one  correspondence  between  lattices  in  the  two  coordinate  systems.  Fig.  5  shows  that  Y  and  W  at  each 

node  [circle  in  Fig.  5]  in  (u,  v)  coordinates  are  directly  mapped  from  Y  and  W ,  respectively,  at  the  corresponding 
node  in  (qd,  qs)  coordinates  without  any  algebraic  manipulation. 

In  Eq.  (26),  a  common  2D  Fourier  argument  U  appears  in  Y  ,  X ,  and  W .  For  each  value  of  U ,  Eq.  (26) 

leads  to  an  over-determined  ID  problem  for  inverse  reconstruction:  7(v)=  fdzW(v,  z)X(z) .  In  order  to 


perform  fast  inversion,  we  invert  the  normal  form  of  the  forward  model:  Y  WT  ~  [WTW  ]  X  for  each  U ,  where 

[  WTW  ]  is  a  MxM  matrix,  with  M the  number  of  layers  in  z  direction.  The  original  W  in  Eq.  (23)  is  a  matrix  with  a 
large  dimension.  The  inverse  problem  now  is  simplified  to  invert  many  (number  of  discrete  value  of  u)  matrices, 
each  with  a  small  dimension  M  The  latter  problem  is  much  more  computationally  efficient  compared  to  the  original 

problem  of  Eq.  (23).  Once  X(u3z)  are  obtained  for  all  ll,  a  2D  inverse  Fourier  transform  produces  Jf(r, z) , 
which  is  the  3D  image  of  optical  parameters  of  the  medium.  Markel  and  Schotland  use  different  procedures  for 
inversion.  In  [13]  a  Fourier-Laplace  inversion  is  applied,  hence,  an  analytic  continuation  of  measured  data  to  the 
complex  plane  is  required  for  the  inverse  Laplace  transform.  In  [14]  an  inverse  procedure  is  performed  in  an 
argument  space,  similar  to  variables  V  here.  Since  V  include  2D  variables,  inversion  in  v  space  could  take  longer 
time  than  that  of  inversion  in  z  space. 
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As  discussed  before,  matrices  W  and  [WTW]  for  each  u  can  be  calculated  in  advance  for  a  uniform 
background  slab  medium.  Assuming  that  a  group  of  experimental  data  has  been  acquired,  die  following  steps  are 
taken  to  produce  a  3D  image  of  the  medium: 

(1)  Obtain  “change”  of  intensities,  7(rd ,  ?s ,  zd ,  z, )  ,  by  subtracting  the  intensity  for  a  uniform  background 
medium  from  the  measured  intensity; 

(2)  Extend  the  (x,  y)  area  and  padding  zeros,  to  overcome  the  wraparound  problem  in  discrete  convolutions; 
[19] 

(3)  Perform  a  dual  2D  fast  Fourier  transform  (FFT)  of  7(rd  ,?s  9zd,zs)  in  the  extended  area  to 
produce  Y(qd,qs9zd,zs)'9 

(4)  Determine  Y  (u,  V,  zd ,  zs )  for  each  u ,  using  a  mapping  procedure  explained  in  Fig.  5; 

(5)  Invert  Y  WT  =  [WTW  ]X  for  each  u ,  which  is  an  inverse  problem  involving  a  MxM  matrix,  with  M 
the  number  of  layers  along  z  direction.  Proper  regularization  according  to  noise  level  needs  to  be  taken  into 
account.  Regularization  will  be  discussed  later  in  the  paper;  and 

(6)  Finally,  perform  an  inverse  2D  FFT  on  X (3,  z)  to  produce  X (r ,  z)  . 

Our  computational  experiments  show  it  takes  only  1-2  minutes  on  a  personal  computer  to  perform  an  inverse 
reconstruction  of  a  3D  image  of  a  medium  with  a  large  number  of  voxels  (for  example,  32x32x20  voxels)  using  this 
HDF  algorithm. 

To  demonstrate  our  concept  of  HDF  tomography  in  3D  image  reconstruction,  an  example  using  simulated  CW 
data  is  presented.  A  slab  turbid  medium,  with  a  transport  mean  free  path  ltr  =  1  mm,  absorption  length  la  =  300  mm, 
and  thickness  zd  =  40  mm,  is  divided  into  20  layers.  A  CW  light  source,  injected  perpendicular  to  the  z*  =  0  plane, 
scans  by  a  2D  32X32  array  on  the  plane,  with  each  pixel  3mmx3mm.  A  2D  array  of  detectors  with  the  same  spacing 
is  located  at  zd  plane  (transmission  geometry).  The  medium,  is  divided  into  32x32x20  voxels,  each  of  dimension 
3x3x2  mm3.  Two  absorbing  objects  are  located  in  the  medium,  each  with  a  volume  3x3x2  mm3.  The  first  one 
located  at  (10,  10,  10)  has  an  absorption  difference  of  A\xa  =  0.01  mm-1  with  the  background.  The  second  one  is 
located  at  (20,  20,  15)  with  an  absorption  difference  of  0.007  mm-1.  The  simulated  data  with  noise  level  of  5%  are 
obtained  using  the  CFM.  The  tomographic  images  are  shown  in  Fig.  6.  As  shown,  the  central  positions  of  3D  image 
of  the  objects  are  correct,  located  at  a  voxel  (10,  10,  10)  with  red  color,  and  a  voxel  (20,  20,  15)  with  yellow  color. 
The  resolution  of  image  is  about  ~  6  mm  in  the  transverse  (x,  y)  plane  and  -  10  mm  along  z  direction.  In  general,  the 
axial  resolution  (along  z  direction)  is  poorer  than  the  lateral  resolution  [on  the  (x,  y)  plane].  In  transmission 
geometry,  two  Green’s  functions  in  the  weight  function  compensate  each  other  when  the  z  position  of  the  object 
changes,  that  leads  to  a  poor  sensitivity  of  the  measured  photon  intensity  to  the  z  position  of  the  object.  The  shapes 
of  3D  image  of  two  objects  are  ellipsoids  with  longer  axis  along  z  direction.  The  absorption  difference  has  the 
maximum  value  at  the  center  of  ellipsoid,  and  decays  gradually  with  increase  distance  from  the  center. 

A  cut-off  in  discrete  lattices  of  qs  and  qd  naturally  introduces  a  kind  of  regularization.  This  regularization  is 
very  effective.  Initial  tests  show  that  even  adding  30%  of  fluctuations  on  simulated  data  of  yr(fd,f5,zrf,z  ) ,  an 
image  similar  to  that  shown  in  Fig.  6  is  still  reconstructed.  The  reason  for  this  is  that  noises  come  from  fluctuations 
at  different  source  and  detector  positions,  which  are  mainly  the  high  frequency  components  of  qs  andqd  .  A  cut-off 

inqs  andqd  naturally  eliminates  these  high  frequency  noises,  such  that  a  stable  image,  especially  in  (x,  y)  plane,  can 
be  reconstructed  in  a  strong  noise  level. 

However,  the  inverse  problem  is  still  ill-posed,  because  contribution  to  the  change  of  intensity  from  a  small 
voxel  deeply  inside  medium  is  weak,  and  is  not  sensitive  to  its  z  position  in  transmission  case.  A  regularization 

procedure  on  inversion  of  Y  WT  =  [WTW]Xi s  still  needed.  The  standard  Tikhonov  regularization  approach  [20] 
is  applied  and  L-curve  [21 ,22]  method  is  used  for  determining  the  best  regularization  parameters. 

This  fast  inverse  algorithm  produces  a  3D  image  in  a  linear  image  regime.  For  nonlinear  image  reconstruction 
procedure,  the  reconstructed  3D  image  provides  a  good  initial  profile  for  further  refining  the  3D  image  taking  the 
nonlinear  effects  into  consideration. 
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The  HDF  inversion  method  can  be  extended  to  a  cylindrical  geometry,  with  an  arbitrary  shape  of  the  (jc,  y) 
cross  section,  for  3D  image  reconstruction.  In  this  geometry,  an  algorithm  using  a  single  Fourier  inversion  has  been 
developed.[23]  This  algorithm  is  limited  to  the  case  that  the  sources  and  the  detectors  are  located  on  a  plane  with 
same  z  coordinates.  The  hybrid-dual-Fourier  inverse  approach  in  cylindrical  geometry  removes  this  restriction,  so 
more  data  can  be  acquired  for  3D  tomography.  The  linear  forward  model  in  cylindrical  geometry  is  given  by 

Y<J&Xs>zd’zs)=  \drdzW(rd,%,r,zd -z,zs -z)X(r,z)  ,  (27) 

where  W (rd ,  rs ,  r ;  zd  -  z,  zs  -  z)  is  the  weight  function,  a  function  of  z d  -  z  and  zs  -  z  due  to  the  1 D  translation 
invariance  of  the  Green’s  function  in  a  homogeneous  background  medium  in  cylindrical  geometry  (assuming 
infinite  z  length).  We  make  a  dual  ID  (along  z  direction)  Fourier  transform  jdz^ dzs e iqdZd e 'q*z*  on  Eq.  (27)  to 
obtain 

r(qd,q„frf,rf)=  J^z#(qd,qi,f,r^,rj)l(qd  +  qs,r),  (28) 

where  Y ,  X ,  and  W  are  the  Fourier  space  quantities  corresponding  that  in  Eq.  (27). 

The  (ID)  linear  hybrid  coordinate  transforms,  u  =  qd+  qSi  and  v  =  qd-  qs,  for  Eq.  (28)  leads  to: 

Y (u, v,rrf,rs)=  jdfW( u, v, r^,r5; x)X (u, r)  ,  (29) 

where  Y  ,  X  ,  and  W  are,  respectively,  Y ,  X  ,  and  W  as  functions  of  U  and  V  .  For  each  value  of  u,  Eq.  (29) 
is  an  over  determined  2D  problem  for  inverse  reconstruction,  namely,  to  determine  a  2D  unknown  value  of 

X(u,r)  from  known  3D  data  of  7(u,  V,?d,rs)  for  each  u.  This  3D-2D  determination  enhances  the  accuracy  of 

3D  image  compared  to  2D-2D  determination  in  the  single-Fourier  transform  inversion.  After  X(u,  r)  are  obtained 
for  all  u,  a  ID  inverse  Fourier  transform  produces  the  image X(r,z)  . 

5.  Discussion 

As  shown  in  Eqs.  (19)  and  (21),  there  is  no  contribution  from  5  wave  to  the  weight  function  for  a  scattering 
object.  This  result  reflects  a  fact  that  no  scattering  effect  exists  for  an  isotropic  angular  distribution.  In  the  regions 
far  from  sources,  the  weight  function  contributed  from  scattering  objects  is  small  because  there  is  no  contribution 
from  the  dominant  5  wave,  as  shown  in  many  results  based  on  the  diffusion  models. [1-5]  This  non-sensitivity  of 
signals  to  the  scattering  objects  deep  inside  the  medium  should  be  considered  in  optical  tomography.  A  pure 
isotropic  distribution  is  never  achieved,  otherwise,  there  will  be  no  flux  in  any  directions.  In  the  diffusive  model,  a 
small  p  wave,  -  (3  /  4;r)Ds  •  VN ,  exists  which  maintains  the  photons  diffusing  to  the  regions  with  fewer  photons. 
The  factor  —  VjV  represents  this  effect.  However,  this  expression  is  valid  only  in  the  regions  where  the  p  wave  is 
much  smaller  than  s  wave,  (  1/4ti)A,  and  does  not  correctly  describe  the  early  photon  propagation  near  sources.  Since 
only  the  weight  function  for  scattering  objects  close  to  sources  plays  an  important  role,  but  it  was  estimated  using 
the  formula  valid  in  regions  far  from  sources,  substantial  error  introduced  in  the  diffusion  forward  model  for 
scattering  objects  is  crucial. 

For  the  weight  function  of  absorbing  objects,  contributions  from  all  spherical  components,  including  5  wave,  are 
given  in  Eqs.  (18)  and  (20).  In  commonly  used  diffusion  formula  the  contribution  from  p  wave  was  neglected.  The 

diffusion  coefficient  originally  derived  in  the  DA  is  D  =  1  /(3//’  +  jua ) ,  that  leads  to 

AD  -  -D(0)“  (3A/JS  *F  A fia  )  .  The  contribution  from p  wave  to  the  weight  function  for  absorbing  objects,  hence, 

should  exist.  But  in  the  later  diffusion  models  A D  is  assigned  only  for  scattering  objects  and  only  s  wave  for 
absorbing  objects  is  taken.  Eqs.  (18)  and  (20)  provide  a  quantitative  estimation  of  weight  function  for  absorbing 
objects  in  regions  close  to  the  source,  as  well  as  far  from  the  source.. 

The  CFM  and  the  HDF  inverse  algorithm  need  further  improvements  in  the  following  aspects.  Further 
improvement  should  be  considered  without  significantly  increasing  complexity  in  computation.  First,  the  second 
cumulant  solution  is  not  accurate  in  the  detailed  shape  of  the  distribution,  especially,  the  front  edge  in  the  Gaussian 
distribution  violates  causality.  An  empirical  distribution,  which  keeps  the  exact  value  of  the  first  and  second 
cumulants,  while  satisfies  the  causality,  can  be  designed  to  replace  the  Gaussian  distribution. 

Second,  the  boundary  condition  is  approximate.  When  a  more  accurate  distribution  /( r,  s,  t)  at  early  time  is 
needed,  the  boundary  condition  for  a  semi-infinite  geometry  should  be 
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I(x,y,z  =  O;0,<f>,t)  =  O  , 


if  cos#  >  0. 


(30) 


This  type  of  the  boundary  condition  was  studied  by  Domke  [24]  for  the  steady  state  case.  The  solution  is  represented 
as  a  superposition  of  a  solution  describing  a  transport  problem  in  an  infinite  medium,  and  a  Fredholm  integral  term, 
which  corrects  this  solution  for  the  appropriate  half-space  boundary  condition.  This  approach  may  be  used  for 
farther  development  of  the  boundaiy  problem. 

Third,  to  consider  the  nonlinear  effects,  /°^’s  in  Eq.  (16)  should  be  replaced  by  the  Green’s  function  in  a  real 
heterogeneous  medium.  Among  the  high-order  perturbative  corrections  of  the  Green’s  function,  the  “self-energy” 
diagram,  which  counts  photon  round  trips  through  a  position  up  to  infinite  times,  plays  an  important  role. 
Gandjbakhche  et  ai  [25]  studied  this  effect  using  a  random  walk  model.  We  find  that  a  renormalization  procedure  for 
this  nonlinear  effect  can  be  performed  after  image  is  obtained  using  a  linear  inversion  process.  This  renormalization 
procedure  can  recover  the  optimal  value  of  the  optical  parameters  and  can  improve  the  resolution  of  image.  The 
detailed  results  of  the  renormalization  will  be  published  elsewhere. 

The  translation  invariance  is  valid  for  the  parallel  geometry  assuming  that  the  (x,  y)  area  is  infinite.  We 
suppose  that  this  assumption  of  the  infinite  area  is  reasonable.  How  much  error  arises  due  to  the  finite  area  of  a 
sample  will  be  studied  in  details. 

Use  of  the  simulated  data  mainly  tests  the  validity  of  the  inverse  algorithm,  does  not  test  accuracy  of  the 
forward  model.  Experimental  data  from  phantoms  and  in  vivo  measurements  in  human  body  will  be  performed  for 
farther  testing  of  our  approach. 

In  summary,  we  have  developed  a  linear  forward  model  of  light  propagation  in  a  turbid  medium  based  on  an 
analytical  cumulant  solution  of  the  radiative  transfer  equation  for  3D  optical  tomography.  The  model  can  be  used  for 
CW,  frequency-domain,  and  time-resolved  measurements  in  parallel  geometries.  This  forward  model  is  more 
accurate  than  the  forward  model  based  on  the  diffusion  approximation  of  RTE.  An  inverse  algorithm  is  developed, 
based  on  a  fast  3D  hybrid-dual-Fourier  tomographic  approach  using  multiple  detectors  and  multiple  sources  in 
parallel  geometries.  This  inverse  algorithm  is  computationally  efficient  and  is  suitable  for  clinical  applications,  such 
as  beast  cancer  detection. 
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Figure  Captions 

Fig.  1  Distribution  function  /( r,  s,  t)  in  an  infinite  uniform  scattering  medium  as  a  function  of  time  t,  using  Henyey- 
Greenstein  phase  function  with  g  =  0.9.  The  detector  is  located  at  R  =  6/fr  =  60  ls  from  the  source  front  along 
direction  of  incident  light,  and  the  direction  is  along  the  incident  direction.  The  solid  curve  is  computed  from 
approximation  up  to  10th  order  of  cumulant;  the  dotted  curve  is  computed  from  approximation  up  to  the  second 
order  of  cumulant,  the  discrete  red  dots  are  from  the  Monte  Carlo  simulation;  the  curve  of  thick  dots  is  from  the 
diffusion  approximation  (DA),  N( r,  t)/4n. 

Fig.  2  The  light  distribution  in  an  infinite  uniform  medium  as  a  function  of  time  at  different  received  angle,  using 
second  cumulant  solution  of  radiative  transfer  equation,  where  detector  is  located,  separately,  at  10  mm  (Fig.  2a)  and 
30  mm  (Fig.  2b)  from  the  source  in  the  incident  direction.  The  parameters  for  this  calculation  are:  1^  =  2  mm,  la  = 
300  mm,  the  phase  function  is  computed  using  Mie  theory  for  polystyrene  spheres  with  diameter  d  =  1.11  pm  in 
water  and  the  wavelength  of  laser  source  X  =  625  nm,  which  gives  the  g-factor  g  =  0.926. 

Fig.  3  A  schematic  diagram  shows  how  to  extend  the  cumulant  solution  of  RTE  from  an  infinite  medium  to  a  semi¬ 
infinite  medium. 

Fig.  4  Backscattered  photon  distribution  /( r,  s  =  -z,  t)  emerging  from  plane  surface  of  a  semi-infinite  turbid 
medium,  as  a  function  of  time,  with  the  source-detector  distance  1  Ilr  on  the  surface  z  =  0  plane.  The  pulse  source  is 
located  at  z  =  0,  incident  along  z  direction.  The  extrapolated  boundary  condition  is  used.  The  solid  curve  is  obtained 
from  cumulant  approximation  (CA),  up  to  the  second  cumulant.  The  dashed  curve  is  from  diffusion  approximation 
(DA).  The  cross  points  are  obtained  from  Monte  Carlo  simulation  (MC). 

Fig.  5  An  example  of  a  6x6  lattice  for  explaining  the  linear  hybrid  transform  from  (qd,  qs)  coordinates  to  (u,  v) 
coordinates. 

Fig.  6.  A  3D  image  reconstructed  using  hybrid  dual  Fourier  tomography.  Two  absorbing  objects,  each  with  the 
volume  3x3x2  mm3,  are  located  inside  a  turbid  medium  with  volume  96x96x40  mm3  divided  into  32x32x20  voxels. 
The  first  one  is  located  at  position  labeled  (10,  10,  10)  with  absoiption  difference  Ap„=  0.01  mm"1.  The  second  one 
is  located  at  position  labeled  (20,  20,  15)  with  absorption  difference  A \i  =  0.007  mm"1.  A  CW  light  source  incident 
perpendicular  to  the  zs  =  0  plane  is  scanned  through  a  2D  32X32  array  at  the  plane,  with  each  pixel  3mmx3mm.  A 
same  sized  2D  array  of  detectors  is  located  at  zd  plane  (transmission  geometry).  The  simulated  data  are  produced 
with  noise  5%.  A  linear  scale  of  color  bar  from  the  maximum  value  to  minimum  value  of  ApCT  is  used.  The  numbers 
labels  the  z  layers  counting  form  source  to  the  detector,  layers  are  separated  by  2  mm. 
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The  anomalous-diffraction  theory  (ADT)  of  extinction  of  light  by  soft  particles  is  shown  to  be  determined  by 
a  statistical  distribution  of  the  geometrical  paths  of  individual  rays  inside  the  particles.  Light  extinction 
depends  on  the  mean  and  the  mean-squared  geometrical  paths  of  the  rays.  Analytical  formulas  for  optical 
efficiencies  from  a  Gaussian  distribution  of  the  geometrical  paths  of  rays  are  derived.  This  Gaussian  ray  ap¬ 
proximation  reduces  to  the  exact  ADT  in  the  intermediate  case  of  light  scattering  for  an  arbitrary  soft  particle 
and  describes  well  the  extinction  of  light  from  a  system  of  randomly  oriented  and  (or)  polydisperse  particles. 
The  implications  for  probing  of  the  sizes  and  shapes  of  particles  by  light  extinction  are  discussed.  ©  2003 
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Anomalous-diffraction  theory  (ADT)  which  was  in¬ 
troduced  by  van  de  Hulst1  for  light  extinction  and 
scattering,  is  one  of  the  simplest  and  most  powerful 
approximations  of  the  interaction  of  electromag¬ 
netic  radiation  with  spherical  and  nonspherical  soft 
particles.  This  approach  has  been  used  in  remote 
sensing  of  cirrus  clouds  and  climate  research,  in 
biophysical  and  biomedical  research,  and  in  other 
applications.2  The  anomalous-diffraction  theory  is 
based  on  the  premise  that  the  extinction  of  light  by 
a  particle  is  primarily  a  result  of  the  interference 
between  the  rays  that  pass  through  the  particle 
with  those  that  do  not.3  This  approximation  is  most 
applicable  to  so-called  soft  particles  with  the  complex 
relative  refractive  index  m  near  1  (\m  -  1|  <5C  1)  and 
with  a  characteristic  dimension  of  size  r  exceeding 
wavelength  A  of  the  incident  radiation  (27 rr/A  >  1)  to 
achieve  a  high  degree  of  accuracy.3-6  This  accuracy 
has  been  observed  to  improve  with  softness  and  non¬ 
sphericity,6  and  with  polydispersity  of  the  particle. 3 

In  this  Letter  we  show  that  ADT  has  a  statistical  in¬ 
terpretation.  The  extinction  of  light  by  particles  mea¬ 
sures  a  probability  distribution  of  the  geometrical  path 
of  the  individual  rays  inside  the  particles  rather  than 
the  sizes  and  shapes  of  individual  particles. 

In  the  framework  of  ADT,1  the  extinction,  ab¬ 
sorption,  and  scattering  efficiencies  of  a  particle  are 
given  by 

Qext  “  -~9t  f(  {1  -  exp [-ikl{mr  -  1)] 
r  p 

X  exp(-klmi)}dP, 

Qabs  =  4:  fl  [1  -  exp(-2klmi)]dP, 

r  p 

Qsca  =  Qext  ~~  Qabs  >  (1) 

where  ?R  represents  the  real  part,  the  wave  number  is 
k  =  2 7?-/ A  for  wavelength  A,  the  complex  relative  re¬ 
fractive  index  is  m  —  mr  -  inti ,  l  is  the  geometrical 


path  of  an  individual  ray  inside  the  particle,  and  P  is 
the  projected  area  of  the  particle  in  the  plane  perpen¬ 
dicular  to  the  incident  light  over  which  the  integration 
is  performed.  The  optical  efficiencies  for  a  system  of 
randomly  oriented  and  (or)  polydisperse  particles  are 
averaged  over  all  the  sizes  and  orientations  of  particles 
weighted  by  their  projection  areas,  i.e., 


Q  * 


Ypq 
Yp  ' 


(2) 


The  integration  in  Eq.  (1)  over  the  projected  area 
for  a  single  particle  at  a  fixed  orientation  or  the  av¬ 
eraging  in  Eq.  (2)  over  the  combined  projected  area 
from  all  sizes  and  orientations  of  particles  can  be  rein¬ 
terpreted  as  an  averaging  over  a  distribution  of  the 
geometrical  path  l  of  rays.  By  dividing  the  (combined) 
projection  area  into  equal-area  elements  and  counting 
the  resultant  geometrical  paths  that  correspond  to  each 
projection  area  element  according  to  their  lengths,  one 
can  find  a  probability  function  p{l)dl  that  describes 
the  probability  that  geometrical  path  l  from  a  ray  is 
within  (l,  l  4-  dZ).  The  probability  function  is  normal¬ 
ized  to  f  p(l)dl  =  1.  By  this  interpretation,  we  can 
rewrite  the  optical  efficiencies  in  Eq.  (1)  as  expected 
values  in  accordance  with  probability  distribution  p(l) 
of  the  geometrical  paths  of  rays.  The  extinction  and 
absorption  efficiencies  in  Eq.  (1)  can  be  expressed  as 

Qext  =  2m  f  {1  -  exp[-ikl(mr  -  1)] 


X  exp (—klrrii)}p(l)  dl , 

Qab*  =  f  [1  -  exp{-2klml)]p{l)dl .  (3) 

Assume  that  the  geometrical  path  distribution  of  rays 
(in  short,  the  ray  distribution)  for  one  particle  with  a 
unit  size  is  p$(l);  then  the  ray  distribution  for  a  par¬ 
ticle  with  the  same  shape,  orientation,  and  a  different 
size  L  is  given  by  p(l)  —  (1  /L)p0(l/L)  from  scaling  of 
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length.  Thus,  a  system  of  such  particles  whose  size  is 
distributed  according  to  a  probability-density  function 
n(x)  has  a  ray  distribution  function 

f  ( l/x)po(l/x)n(x)x2dx 
Ppol(Z)  ~  7  ,  (4) 

/  n{x)x2  dx 

weighted  by  the  projection  area  of  individual  particles 
that  is  proportional  to  x2.  The  subscript  pol  or  m  is 
used  to  denote  a  polydispersed  particle  or  one  that  is 
randomly  oriented,  respectively. 

Let  us  consider  a  unit  spheroid  with  a  semisize  6  =  1 
of  the  revolutional  axis  and  an  axial  ratio  e  and  with 
an  angle  x  between  the  propagation  direction  of  the 
incident  beam  and  the  revolutional  axis  of  the  spher¬ 
oid.  The  geometrical  length  of  a  ray  and  the  projection 
area  for  such  a  spheroid  have  been  calculated.5  The 
geometrical  path  distribution  of  the  rays  can  then  be 
found: 

Po(l )  =  \  (e~2  sin2 x  +  cos2  x)W{l) 


x  gr  2 

(e  2  sin2  x  +  cos2  x )1/2 


(5) 


where  H(x)  is  a  Heaviside  function.  The  ray  distribu¬ 
tion  for  a  system  of  such  spheroids  at  a  fixed  orienta¬ 
tion  x  with  a  log-normal  size  distribution,7 


n(x)  = 


(27r)1/2cr  r 


exp 


f  ln2(r  /  am)~\ 
2a-2 


(6) 


is  given  by 

_  /n  (f  "2  sin2  x  +  cos2  x)l 

Ppol(Z)  . 


erfc{(l/V2o-)ln[(6  2  sin2  x  +  cos2  ;y)1/2Z/2am]} 

am2  exp(2cr2)  1  ‘ 

from  Eq.  (4),  where  erfc(x)  is  the  complementary  error 
function.  The  ray  distribution  becomes 


out  by  the  averaging  over  the  polydispersity  and  the 
orientation  of  the  particle.  The  shape  characteristics 
of  an  individual  particle  are  expected  to  be  further 
washed  out  if  particles  of  different  shapes  are  involved. 
Thus  the  ray  distribution  p(l)  of  a  system  of  particles 
such  as  a  bacterial  suspension,  biological  cells,  or  cir¬ 
rus  clouds  where  particles  are  polydisperse,  randomly 
oriented,  and  (or)  of  multiple  shapes  approaches 
a  probability-density  function  p(l)  that  is  charac¬ 
terized  essentially  by  the  mean  geometrical  path 
(l)  —  [  lp(l)dl  and  the  mean-squared  geometrical  path 
(l2)  =  [  l2p(l)dl  of  rays  inside  the  particles.  One 
natural  choice  of  p(l)  here  is  the  Gaussian  probability- 
distribution  function,  which  follows  the  same  spirit 
as  the  well-known  central -limit  theorem. 8  We  should 
point  out  that  this  choice  does  not  satisfy  p(l  <  0)  —  0, 
but  the  contribution  from  near  the  Z  =  0  region  in  the 
ray  distribution  is  much  smaller  than  that  from  other 
regions  and  hence  can  be  ignored. 

Let  us  now  assume  that  the  ray  distribution  is  given 
by  a  Gaussian  distribution: 


p(x)  = 


(*-  m)21 

2a-2 


(9) 


The  extinction  and  scattering  efficiencies  are  then 
given  by 


Qext  =  2-2  cos[&(??2r  1)  (//,  —  kcr2rrii)] 

„  f  7  k2cr2[(mr  ~  l)2  -  rrii 2]] 

x  exp '-kpmi - - - - — 'i 


Qabs  —  1  —  exp[-2  klUiip  —  klYli  O'2)] 


GO) 


from  Eqs.  (3)  after  a  straightforward  integration. 
The  optical  efficiencies  [Eqs.  (10)],  in  the  intermediate 
case  limit  [k(mr  -  1)1  «  1  and  krriil  «  1,  where  Z  is 
the  geometrical  path],9  reduce  to 


Qext  =  2kn i,il)  +  k\(mr  -  l)2  -  m,-2](Z2), 
Qabs  =  2krrti{l)  —  2k2m2{l2) , 


Ppol,  rn  (Z) 

f  pVQ\(l)ir e2{e~2  sin2  x  +  cos2  ^)1/2d  cos  x 


I  7 re2(e  2  sin2  x  +  c°s2  ^)1//2d  cos  x 
o 

for  such  particles  randomly  oriented  where  the 
projection  area  of  the  cylinder  is  proportional  to 
7 T€2(e~2  sin2  x  +  cos2  x )1/2-5  It  is  worth  noting  here 
that  the  ray  distribution  for  a  simple  spheroid  at  a 
fixed  orientation  [Eq.  (5)]  is  triangular,  regardless  of 
the  axial  ratio  of  the  spheroid.  This  fundamental  geo¬ 
metrical  characteristics  facilitates  a  simple  rescaling 
of  the  radius  to  calculate  the  optical  efficiencies  from 
a  sphere  for  a  spheroid.5 

The  ray  distributions  from  a  single  spheroid,  a  single 
randomly  oriented  spheroid,  a  system  of  polydisperse 
spheroids  at  a  fixed  orientation,  and  a  system  of 
randomly  oriented  polydisperse  spheroids  are  plotted 
in  Fig.  1.  It  is  clear  from  the  figure  that  the  shape 
characteristics  of  an  individual  particle  are  washed 


Qsca  —  k2\ m  -  1|2(Z2),  (11) 

where  the  mean  and  the  mean -squared  geometri¬ 
cal  paths  are  given  by  <Z)  =  p  and  <Z2}  =  p2  +  a2, 


Ray  Path  Ray  Path 

(a)  (b) 


Fig.  1.  Ray  distributions  for  a  spheroid  at  a  fixed  orien¬ 
tation  x  =  0  (FX),  randomly  oriented  (RN),  polydisperse 
at  a  fixed  orientation  (POL  FX),  and  randomly  oriented 
polydisperse  (POL  RN).  The  axial  ratio  of  the  spheroid  is 
(a)  e  =  2  and  (b)  e  =  0.5.  Log-normal  size  distribution 
rc(x)  with  am  —  1  and  a  =  0.2  is  also  plotted  as  insets. 
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(a)  (b) 


Fig.  2.  Extinction  and  absorption  efficiencies  of  (a)  a 
sphere  and  (b)  a  polydisperse  sphere  with  a  log-normal 
radius  distribution  of  am  =  1  and  a  =  0.2  calculated  with 
Mie,  ADT,  and  Gaussian  ray  approximations.  Complex 
refractive  index,  m  —  1.05  -  i0.0005.  The  size  distribu¬ 
tion  has  already  been  shown  as  insets  in  Fig.  1. 

respectively.  These  results  agree  exactly  with  those 
for  the  intermediate  region  over  which  the  Rayleigh  - 
Gans  approximation  and  the  anomalous-diffraction 
approximation  of  light  scattering  from  small  particles 
overlap.1*9  This  means  that  Eqs.  (10)  from  our  Gauss¬ 
ian  ray  approximation  reduce  to  the  exact  ADT  in  the 
intermediate  case. 

Figure  2  compares  the  extinction  and  absorption  ef¬ 
ficiencies  calculated  by  the  exact  Mie  theory,  the  exact 
ADT  [Eqs.  (1)  and  (3)],  and  our  Gaussian  ray  approxi¬ 
mation  [Eqs.  (10)]  for  a  weakly  absorbing  sphere  and 
a  system  of  the  same  spheres  with  a  log-normal  radius 
distribution  [Eq.  (6)]  of  am  =  1  and  a  =  0.2.  Both 
our  Gaussian  ray  approximation  and  the  ADT,  unlike 
the  exact  Mie  calculation,  tend  to  underestimate  the 
optical  efficiencies.  This  fact  is  well  known. 10,11  The 
absorption  efficiency  from  our  Gaussian  ray  approxi¬ 
mation  agrees  extremely  well  with  the  ADT;  at  most 
it  differs  by  2%  from  the  exact  Mie  calculation  in  this 
comparison.  The  extinction  efficiency  agrees  well 
with  the  exact  Mie  calculation  in  the  intermediate  re¬ 
gion  for  both  single  spheres  and  polydisperse  spheres, 
as  expected.  The  Gaussian  ray  approximation  for 
the  polydisperse  spheres  approaches  the  exact  ADT 
calculation  with  maximum  relative  errors  of  3.5% 
compared  to  the  ADT  and  of  7%  compared  to  Mie 
theory. 

From  our  statistical  analysis  of  the  anomalous- 
diffraction  theory  of  light  extinction,  light  extinction 
depends  solely  on  the  probability  distribution  of  the 
geometrical  paths  of  individual  rays  inside  the 
particles  rather  than  on  the  size  or  the  shape  of  an  in¬ 
dividual  particle.  Thus  the  optical  efficiency  equiva¬ 
lence12  can  easily  be  achieved  from  different-shaped 
particles  or  particles  of  different  size  distributions 
as  long  as  they  share  a  common  geometrical  path 
distribution  of  rays. 

The  geometrical  path  distribution  of  rays  can  be 
approximated  by  a  Gaussian  probability  distribution 
function  for  a  system  of  particles  in  which  the  par¬ 


ticles  are  randomly  oriented,  polydisperse,  and  (or) 
multiple  shaped.  For  such  a  system  of  particles  the 
light-extinction  measurements  essentially  determine 
the  mean  and  the  mean-squared  geometrical  paths 
of  rays  from  all  particles  in  the  system.  The  shape 
and  size  of  an  individual  particle  can  be  deduced 
only  with  a  priori  information  on  the  shape  and  (or) 
the  size  distribution  of  the  particles  involved.  The 
pursuit  of  the  mean  and  the  mean-squared  paths 
from  fitting  Eqs.  (10)  to  experimental  data,  or  the 
general  geometrical  path  distribution  of  rays  p{l)  of 
particles  from  solving  the  inverse  problem  in  Eqs.  (3), 
provides  an  alternative  approach  to  particle  sizing 
and  shaping.  We  note  that  we  have  restricted  this 
study  to  extinction  of  light  from  particles  of  the 
same  type  (a  common  refractive  index,  m).  This 
statistical  interpretation  of  ADT  opens  a  new  way 
to  calculate  optical  efficiencies  of  soft  particles  of 
different  shapes  by  use  of  the  probability  distribution 
of  the  geometrical  paths  of  individual  rays  inside 
particles. 
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We  have  developed  a  hybrid-dual-Fourier  tomographic  algorithm  for  a  fast  3D  image  recon¬ 
struction  in  turbid  media  using  multiple  sources  and  multiple  detectors.  This  algorithm  can  be 
applied  to  different  forward  models  based  on  a  diffusion  approximation*  or  a  cumulant  analytical 
solution  of  the  radiative  transfer  equation2’3  in  CW,  frequency-,  and  time-domains  as  long  as  the 
translational  invariance  in  the  lateral  direction  exists.  It  provides  a  general  scheme  of  fast  numerical 
three-dimensional  image  reconstruction  using  multiple  sources  and  detectors  and  is  different  from 
previous  diffraction  tomography  approaches.1,4 

By  scanning  a  point  source  on  the  grids  of  the  input  plane  of  a  slab,  the  measured  light  intensity 
I(pd-Ps)  on  a  detector  array  on  the  exit  plane  of  the  slab  may  form  a  4D  data  set  where  pd  and  ps 
are  the  coordinates  on  the  input  and  exit  plane.  An  additional  dimension  of  the  data  set  may  be 
added  in  frequency-modulated  or  time-resolved  measurements.4,5  Based  on  the  general  properties, 
of  the  translational  invariance  along  the  lateral  direction,  a  dual  two  dimensional  Fourier  transform 
on  both  lateral  coordinates  of  the  detector  pd  and  the  source  ps  is  applied  to  the  intensity  data  set. 
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A  hybrid  transform  is  then  performed  on  the  spatial  frequencies  to  rotate  the  coordinate  system  in 
the  Fourier  space.  This  data  rearrangement  results  in  an  inversion  problem  in  the  form  of: 

A/(u,  v,2d,zs)  =  J  dzW{u,v,z,zd,zs)X(xi,z)  (1) 

in  which  A /  is  the  change  of  intensity  in  the  rotated  Fourier  space  (u,  v),  W  is  the  weight  function 
and  X  is  the  Fourier  quantity  of  the  absorption  and/or  scattering  deviation  from  the  uniform 
background.  The  integral  is  performed  on  the  axial  coordinate  z  of  the  inhomogeneity,  and  zd  and 
zs  are  the  axial  coordinates  of  the  detector  and  the  source. 

The  reconstruction  is  then  performed  for  each  fixed  value  of  the  spatial  frequency  u  in  Eq.  (l). 
A  series  of  such  a  one  dimensional  inversion  over  the  axial  coordinate  z  is  carried  out  to  obtain 
the  Fourier  quantity  X(u,z)  for  each  u.  The  inverse  Fourier  transform  of  X(u,z)  provides  a  three 
dimensional  reconstruction  of  the  inhomogeneities  in  the  turbid  medium. 

The  major  advantage  of  this  approach  is  its  speed  and  tolerance  to  noise.  The  reduction  of  the 
image  reconstruction  to  a  series  of  one  dimensional  inverse  problem  of  the  dimension  of  the  number 
of  discrete  divisions  in  the  axial  direction  from  our  data  rearrangement  very  much  reduces  the 
complexity  of  the  image  reconstruction  problem.  The  number  of  operations  required  increases  only 
linearly  with  the  total  number  of  detector-source  pairs.  A  complete  reconstruction  of  a  32  x  32  x  20 
volume  is  within  few  minutes  using  a  1GHz  CPU.  The  tolerance  to  noise,  in  particular,  the  additive 
noise  is  achieved  due  to  the  different  spatial  frequency  spectrum  of  the  signal  and  the  noise.  It  is 
straightforward  in  the  Fourier  space  to  filter  out  noise  which  usually  appears  in  higher  frequency 
components. 

We  will  present  the  reconstruction  results  for  a  CW  and  a  time-resolved  reconstruction  of  a  slab 
with  a  thickness  of  40mm  and  a  similar  optical  property  of  a  human  breast.  We  will  also  discuss  in- 
detail  the  reconstruction  algorithm,  especially,  the  choice  of  the  regularization  parameters  for  Eq.  (1) 
from  a  statistical  analysis  of  the  reconstruction  procedure.  The  merit  of  applying  this  algorithm  to 
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CW,  time-  and  frequency- domains  for  imaging  will  also  be  discussed. 

This  work  is  supported  by  US  Army  Medical  Command.  One  of  the  author  (M.  Xu) 
the  support  of  the  Department  of  the  Army  (Grant#  DAMD17-02-0516). 
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ABSTRACT 

We  report  on  the  effect  of  the  nonlinear  multiple  passage  on  optical  imaging  of  an  absorption  inhomo¬ 
geneity  of  finite  size  deep  inside  a  turbid  medium  based  on  a  cumulant  solution  to  radiative  transfer. 
An  analytical  expression  for  the  nonlinear  correction  factor  is  derived.  Comparison  to  Monte  Carlo 
simulations  reveals  an  excellent  agreement.  The  implication  on  optical  imaging  is  discussed. 

Keywords:  nonlinear  correction,  multiple  passage,  radiative  transfer,  optical  imaging 

1.  INTRODUCTION 

The  principle  of  optical  imaging  of  turbid  media  (such  as  tissues)  is  to  locate  and  reconstruct  the  optical 
properties  (absorption  and  scattering  coefficients)  of  embedded  inhomogeneities  (such  as  tumor)  in  the 
hope  of  identification  by  inverting  the  difference  in  time-resolved  or  frequency-modulated  photon  trans¬ 
mittance  due  to  the  presence  of  the  inhomogeneities  through  either  iterative  or  noniterative  methods. 
The  key  quantity  involved  is  the  weight  function  which  quantifies  the  influence  on  the  detected  signal 
due  to  the  change  of  the  optical  parameters  of  the  medium.  The  diffusion  approximation  to  radiative 
transfer  provides  an  adequate  model  for  the  weight  function  (or  Jacobian)  for  a  small  and  weak  ab¬ 
sorption  inhomogeneity  far  away  from  both  the  source  and  the  detector.  However,  the  weight  function 
predicted  by  the  linear  perturbation  approaches  is  no  longer  valid  when  the  absorption  strength  is  not 
small.  This  can  be  attributed  to  the  multiple  passage  of  a  photon  through  one  single  abnormal  site. 

The  change  of  the  light  intensity  A7  at  the  detector  rd  due  to  the  presence  of  an  absorption  site  at 
r  from  a  modulated  point  source  at  rs  is  expressed  as 

A/  =  -6/j.aVG(rd,  w|r)G(r,  w|rs)  (1) 

to  the  first  order  of  Born  approximation  where  S/j,a  is  the  excess  absorption  of  the  absorption  site  whose 
volume  is  V ,  d  is  the  modulation  frequency  of  light,  and  G  is  the  propagator  of  photon  migration  in 
the  background  medium.  Here,  the  Green’s  function  G(r2,m|n),  in  general,  depends  on  the  detail  of 
light  scattering  inside  the  medium,  and  the  incident  and  outgoing  directions  of  light. 

When  the  absoiption  strength  is  not  small  (Sfj,aV  ^  1),  photon  loss  due  to  multiple  passage  of  the 
absorption  site  is  appreciable  and  can  not  be  ignored.  The  expression  for  A I  in  Eq.  (1)  needs  to  be 
modified  to  include  the  contributions  from  multiple  visits  of  the  site  by  the  photon.  Fig.  (1)  illustrates 
the  most  important  correction  (a  ‘’self-energy”  correction)  which  takes  into  account  the  repeated  visits 
made  by  a  photon  to  the  site  up  to  an  infinite  times. 

Further  author  information:  (Send  correspondence  to  M.  Xu) 
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Figure  1.  Self- energy  correction  to  the  multiple  passage  effect  on  light  absorption. 

Assuming  that  the  center  of  the  absorption  site  is  located  at  r  and  far  away  from  both  the  source 
and  the  detector,  the  change  of  the  detected  light,  A/,  is  now  given  by 


A/  -  -G(rd,uj\r)V6fj,a(r)  ^  f“iVself(a;;i?)V^/ia(r)lnG(r,a;|rs) 

n~0 

VSfXa(  f) 


(2) 


=  -G(rd,cu|r) 


1  +  -^seif(w;  R)V6/j,a(T) 


crG(r,u;|rs) 


where 

iVseif  (w:  R)  =  ^  Jv  Jv  G{ r2,  u/|ri  )d3r2dzr1  (3) 

is  the  self-propagator  which  describes  the  probability  that  a  photon  revisits  the  volume  V  of  size  R. 
Here  G(r^,o;|r)  and  G(r ,  w|rs)  are  well  modelled  by  the  center-moved  diffusion  model  as  long  as  the 
separations  |rd  -  r| ,  |rs  -  r|  »  lt  where  lt  is  the  transport  mean  free  path  of  light  in  the  medium.2 
However,  the  diffusion  Green’s  function  can  not  be  used  in  Eq.  (3)  to  evaluate  Nseu(uj:R)  where  rx  is 
in  the  proximity  of  r2-  By  comparing  Eq.  (2)  to  Eq.  (1),  the  nonlinear  multiple  passgage  effect  of  an 
absorption  site  can  be  summarized  by  the  nonlinear  correction  factor  [l  +  lVself(u/;  R)VSfia(r)}  This 
factor  serves  as  a  universal  measure  of  the  nonlinear  multiple  passage  effect  as  long  as  the  absorption 
site  is  far  away  from  both  the  source  and  the  detector  and  its  size  is  much  smaller  than  its  distance  to 
both  the  source  and  the  detector. 

In  this  article,  we  will  derive  an  analytical  expression  for  the  self-propagator  to  understand  the 
nonlineai  multiple  passage  effect  on  light  absorption  using  our  cumulant  solution  to  radiative  transfer. 
The  nonlinear  correction  factor  [l  +  A'self  R)V 5 /ua(r)]  1  of  our  result  is  shown  to  be  in  an  excellent 
agreement  with  the  Monte  Carlo  simulations  for  continuous  wave  light. 


2.  THEORY 

To  take  into  account  the  higher  order  contributions  from  the  absorption  inhomogeneity,  the  behavior 
of  the  photon  migration  in  a  short  distance  must  be  considered.  Although  the  photon  distribution  is 
almost  isotropic  at  an  absorption  site  deep  inside  the  medium,  the  diffusion  approximation  is  still  not 
appropriate  here.  The  separation  between  the  two  points  rj  and  r2  within  the  volume  in  Eq.  (3)  is 
small.  The  photon  propagator  Y(r2,t|r1,s),  which  represents  the  probability  that  a  photon  propagates 
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from  position  ri  with  propagation  direction  s  to  position  r2  in  time  t,  when  r2  is  in  the  proximity  of 
ri>  is  governed  by  the  radiative  transfer  equation  rather  than  the  diffusion  equation. 

Recently  we  have  shown  that  the  propagation  of  photon  inside  a  turbid  medium  (the  radiative 
transfer  equation)  can  be  solved  analytically  using  a  cumulant  expansion  of  the  photon  distribution 
function.  The  propagation  of  photon  was  found  to  transform  from  an  initial  ballistic  motion  at  early 
time  and  then  gradually  to  a  center-adjusted  diffusion  at  later  time.  The  propagator  of  photon  density 
(the  Green’s  function)  in  an  infinite  uniform  medium  is  given  by4 


lV(r,t|r0,s0)  = 


[47TD(t)t]3/2 


exp  - 


(r  -  r0  -  s0A(t))2 
4D(t)t 


ignoring  the  small  difference  in  the  diffusion  coefficient  along  different  directions  where  the  absorption 
coefficient  is  pa,  the  time-dependent  diffusion  coefficient  is 


~  zl  ^ ~  t1  “  exp (-ct/k)]  -  ^  [1  -  exp (-ct/lt)}2^  (5) 

and 

A(t)  =  lt[  1  -  exp(—ct/lt)]  (6) 

is  the  average  center  of  photons  which  moves  with  speed  c  initially  and  approaches  the  transport  mean 
free  path  lt  in  the  long  time  limit.  The  Green’s  function  for  parallel  geometries  can  be  obtained  by  the 
method  of  image  sources.0 


2.1.  Propagator  of  an  isotropic  point  source 

Let’s  now  consider  the  propagator  lV(r,i|ro,so)  at  the  inhomogeneity  site  ro  =  0  (the  origin  of  space) 
deep  inside  the  medium.  The  photon  distribution  at  r0  is  almost  isotropic,  but  is  anisotropic  scattering. 
The  effective  propagator  can  then  be  obtained  by  averaging  (4)  over  the  propagation  direction  s0  of 
light  over  the  47r  solid  angle,  and  is  given  by  [see  Appendix  A] 


X  (exp  [-(LZ^))2]  _  exp  \JZ±M 
\  1  4Z?(t)t  P  AD(t)t 


(4-7t)3/2  (D(t)t)  1/2rA  (t) 

exp  ^))2n 

P  [  4 D(t)t  J  J  • 


This  reduces  to 


Aeff(r,  t)  =  ^^.r2 —5{r  ~  ct),  for  t  ->•  0H 


Nea(r  t)  = _ exp(  Map  f  (r  ~  It)2 

’  (47r)3/2(D0C^)i/2r/f  iexP  4Dsot 


(r  +  h)' 
4D-~t 


,  for  t  »  1 


in  early  and  late  time  limits  where  Dx  =  ltc/ 3. 

The  temporal  Fourier  transforms  of  the  asymptotic  equations  (8)  and  (9)  are  given  by 


-Veff(r,w) 


1  lr  Nrl 

4iAexp 
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and 


N* (r»w)  =  ^Brnk  [exp(-K  lr  “  **l)  ~  exp(-K(r  +  lt))]  (11) 

respectively,  where  k  =  i/3 (jTa  -  ico)/ltc  whose  sign  is  chosen  with  a  nonnegative  real  part.  In  the  limit 
of  small  k  <C  1,  Eq.  (11)  simplifies  to 

lit  ■  <12> 

This  is  the  case,  for  example,  that  a  continuous  wave  propagates  in  a  nonabsorbing  medium.  The 
errornous  divergence  at  the  zero  separation  in  the  diffuse  Green’s  function 


exp  {—nr) 
47rD00r 


(13) 


is  removed  in  our  formulation  of  the  propagation  of  an  isotropic  point  source. 

The  asymptotic  equation  (11)  from  the  late  time  limit  provides  a  good  approximation  for  A'etf(r,u>) 
when  r>  lt  [see  Fig.  (2)].  The  contribution  to  JVeff(r\w)  when  r  <  lt  is  from  either  ballistic  or  diffusive 
photons,  hence  an  improvement  to  Eq.  (10)  can  be  made 


^  4^eXp  ^  s^^’  r<lt  (14) 

to  include  the  contribution  from  diffusive  photons.  The  effective  propagator  in  temporal  Fourier  space 
iVeff(r,w  =  0)  and  its  asymptotic  behaviors  (10),  (11)  and  (14)  axe  shown  in  Fig.  (2).  The  diffusion 
Green's  function  has  a  huge  error  for  small  r. 


2.2.  Self  propagator  for  a  finite  volume 

For  an  absorption  site  of  a  finite  volume  V  deep  inside  the  medium,  say  a  sphere  of  radius  R  C  L 
where  L  is  the  dimension  of  the  medium,  the  self-propagator  Nsel{(t;R )  for  this  volume  which  denotes 
a  photon  revisits  the  site  in  time  t  is  written  as: 

Nseii(t‘,  R )  =  y2  fvfv  NeS (N  -  ra  | ,  t)d3rid3r2 

1  f2R 

=  y  I  ATeff(r,f)7o(r)47rr2dr  (15) 

where 


70  (r) 


3r  1  (  r  \  3 
4R  +  16  R 


(16) 


is  the  characteristic  function  for  a  uniform  sphere.6"8  This  characteristic  function  has  a  form  of 


7o(r)  =  1  -  (S/4V>  +  . . . 


(17) 


for  an  arbitrary  particle  where  S  is  the  surface  area  of  the  particle.  This  self  propagator  (15)  for  a 
finite  volume  is  quite  different  from  the  self-propagator  of  a  point,  obtained  by  setting  r  —  0  in  (4)  or 
(7),  i.e., 


A'eff(O.t)  = 


exp  (-Mat) 
(4*Z>(t)t)3/2  exp 


A (f)2  ’ 

4 D{t)t  ’ 


(f  >  0). 
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Figure  2.  The  effective  propagator  in  temporal  Fourier  space  Nea(r,u>  =  0)  for  photon  migration  in  a  nonab- 
sorbmg  medium.  Its  approximations  by  (14)  when  r  <  lt  and  by  Eq.  (11)  when  r  >  lL  are  also  plotted.  The 
diffusion  Green’s  function  has  a  huge  error  for  small  r. 


See  Fig.  (3).  This  difference  comes  from  the  fact  that  Eq.  (15)  includes  the  contribution  from  the 
ballistic  motion  of  the  photon  when  the  photon  flies  across  the  site  while  Eq.  (18)  does  not  contain  this 
effect.  The  former  manifests  itself  in  Fig.  (3a)  as  the  linear  decay  of  NseS(t:  R)V  in  the  form  of  7o(ct) 
near  the  origin. 

The  self-propagator  in  temporal  Fourier  space  is  thus  obtained  by  a  temporal  Fourier  transform  of 


r+oo 

=  J  iVself(t;  R)  exp(iu>t)dt  (19) 

1  f2^  9  r+oo 

=  y  Jo  drj0(r)4nr2 1  dtNeS(r,t)exp(iu;t) 

1  f2R 

=  y  jo  •/Veff(r,o;)7o(r)47rr2dr. 

The  lower  limit  of  integration  is  0+,  emphasizing  that  t  =  0  should  be  excluded  from  integration.  Note 
lim{_0+  Aeff(r.t)  =  0  for  our  cumulant  photon  density  function.  This  is  not  the  case  for  the  diffusion 
Green’s  function.  A  numerical  quadrature  is  generally  required  to  compute  this  self  propagator  (19). 
A  crude  estimation  of  AW  (o';  R)  can  be  obtained  from  the  asymptotic  behavior  (11)  and  (14)  of 

-WWW-  i-e-- 


AseffW  R)  - 


1  rmin(2R.lt)  j  r  r-i 

VJo  ^J^exp  (iu  -  Iia)-' 'n(r)4m*dr  (20) 
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(a)  (b) 


Figure  3.  The  self-propagators  for  a  finite  volume  and  a  point:  (a)  Nsel{(t;  R)  and  (b)  Neii(0,t). 


This  reduces  to 

JWu  =  0;/l)  =  £ 

for  a  continuous  wave  propagating  inside  a  nonabsorbing  medium  (u  =  =  k  =  0).  This  estimation 

turns  out  to  be  amazingly  good.  Fig.  (4)  plots  iVself(u>  =  0;  R)  from  numerical  quadrature  and  the 
crude  estimation  (21). 

3.  RESULTS  AND  DISCUSSION 

The  multiple  passage  effect  due  to  the  absorption  site  can  now  be  computed  using  the  self-propagator 
Eq.  (19)  derived  here.  For  large  sites,  the  self-propagator  Nse\ f(o>  =  0;  R)  increases  inverse  proportional 
to  its  size  (iVseif  oc  R  1)  from  Eq.  (21);  hence  the  nonlinear  correction  factor  has  a  form  of 

- = - - - ~  q  + 

l  +  Nself(uj:R)V8»a(r)  5  ltc  ’  {22) 

dependent  on  the  area  of  the  absorption  site  for  large  R. 

Monte  Carlo  methods  have  been  extensively  used  in  simulation  of  photon  migration.9- 10  We  perform 
Monte  Carlo  simulations  on  a  uniform  nonabsorbing  and  isotropic  scattering  slab  (the  anisotropic  factor 
of  scattering  g  —  0).  The  units  of  length  of  time  are  chosen  such  that  the  mean  scattering  length 
^  “  1/ft  =  1  and  the  speed  of  light  c  =  1.  The  transport  mean  free  path  is  hence  It  =  1  and  the 
thickness  of  the  slab  is  assumed  L  —  80^.  An  absorption  spherical  site  is  located  at  the  center  of 
the  slab  (0.  0.  L/2)  with  radius  R  whose  absorption  and  scattering  coefficients  are  /ia.2  —  8(ia  =  0.01 


M 

4c 


rt3 


R  <  h/ 2 


384i?5+160i2fi3-60i3fl2+3/?  „  _  , 

- 320 mrc  -  R>lt/2 
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Figure  4.  The  self  propagator  Ngcui-^  —  0;  R)  and  its  estimator.  The  diffusion  self-propagator  for  continuous 
waves  is  also  plotted. 


and  fxa> 2  =  fis  respectively.  The  photon  is  incident  at  the  origin  on  the  left  boundary  of  the  slab 
z  =  0  in  the  normal  direction  of  the  surface.  Each  photon  is  traced  until  it  escapes  the  slab  through 
either  the  left  or  the  right  boundary.  The  correlated  sampling  is  used  in  simulation  to  reduce  variance. 
A  single  simulation  is  used  to  compute  the  emitted  photon  density  70  for  the  uniform  background 
(nonabsorption  slab)  and  I  for  the  slab  with  the  absorption  site  present. 

The  nonlinear  correction  factor  fl  +  Nse\{(oj',  R)V <5/ia(r)]  1  in  Eq.  (22)  can  be  extracted  from  the 
change  of  the  detected  light  intensity  due  to  the  presence  of  the  absorption  site  in  Monte  Carlo  sim¬ 
ulations  according  to  Eq.  (2).  Fig.  (5)  plots  the  theoretical  nonlinear  correction  factor  and  that  from 
Monte  Carlo  simulations.  “Back”  and  “Forward”  denote  the  cases  where  light  emits  from  the  left 
(z— 0)  and  the  right  (z  =  L)  boundaries,  respectively.  The  agreement  between  our  theoretical  result 
and  Monte  Carlo  simulations  is  excellent  except  for  extremely  small  sizes  of  inhomogeneities. 

Figs.  (6)  and  (7)  plot  the  nonlinear  correction  factor  versus  the  variation  of  the  modulation  frequency 
of  light  for  a  fixed  absorption  strength  and  versus  the  size  of  the  absorption  site  with  a  fixed  modulation 
frequency  of  fight  respectively.  With  the  increase  of  the  modulation  frequency  of  fight,  the  nonlinear 
correction  becomes  less  accentuated.  The  dependence  on  the  size  of  the  inhomogeneity  is  no  longer 
monotonic  for  modulated  fight  while  the  nonlinear  correction  factor  decreases  monotonically  with  the 
increase  of  the  size  for  continuous  wave  fight.  The  phase  delay  is  in  the  order  of  a  few  degrees  in  the 
cases  investigated. 

The  typical  value  of  the  absorption  coefficient  of  human  tissues  is  around  O.OOlps""1  while  the 
scattering  coefficient  is  about  lps'1.  Hence  the  absorption  and  scattering  ratio  is  in  the  order  of  0.001. 
This  should  be  compared  to  our  results  listed  here  where  the  corresponding  ratio  is  0.01  and  one  order 
of  magnitude  stronger.  The  nonlinear  correction  factor  for  absorption  inhomogeneities  such  as  tumors 
in  human  tissues  is  not  appreciable  unless  the  size  of  the  inhomogeneity  is  R  ~  5 It  or  larger. 
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Figure  5.  The  nonlinear  correction  factor  from  the  theoretical  self-propagator  Eq.  (19)  and  Monte  Carlo  simu¬ 
lations.  “Back”  and  “Forward”  denote  light  emitting  from  the  left  (z=0)and  the  right  (z  =  L)  boundaries.  The 
excess  absorption  is  5(ia  =  0.01. 


Figure  6.  The  nonlinear  correction  factor  versus  the  variation  of  the  modulation  frequency  of  light.  The  size  of 
the  absorption  sphere  is  R  =  3 lt.  The  excess  absorption  is  djj.fJ  =  0.01. 


R/l, 


Figure  7.  The  nonlinear  correction  factor  versus  the  variation  of  the  size  of  the  absorption  site.  The  modulation 
frequency  of  light  is  u  -  0.1.  The  excess  absorption  is  Sfia  =  0.01. 

In  conclusion,  we  have  derived  an  analytical  expression  for  the  nonlinear  correction  factor  which 
agrees  well  with  Monte  Carlo  simulations.  The  effect  of  the  nonlinear  multiple  passage  of  an  absorption 
site  on  optical  imaging  only  becomes  appreciable  when  the  size  of  the  inhomogeneity  is  5 lt  or  larger 
for  human  tissues. 


APPENDIX  A.  DERIVATION  OF  NEFF(R,T ) 

The  spatial  Fourier  transform  of  (4)  is  given  by 


AT(k,f|r0,so)  =  J d3r exp(— ik  •  r)iV(r,  f |ro,  s0)  =  exp  ( -k2D(t)t  —  (j,at  -  fk  •  s0A(f)^  .  (23) 

Hence,  the  effective  propagator  in  spatial  Fourier  space  at  r0  is  expressed  as 

•Neff (k, t)  =  /  d2soN (k,  f |r0, s0)  =  exp  (~k2D{t)t  -  fxat\  (24) 

by  averaging  (23)  over  the  propagation  direction  s0  of  light  over  the  4x  solid  angle.  The  effective 
propagator  in  real  space  is  then  obtained  by  an  inverse  spatial  Fourier  transform  0f  (24): 


-Neff  (r,  t)  =  /  exP(zk  •  r)exP  l-k2D(t)t  - 
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1.  Introduction 


Over  the  past  decade,  the  study  of  optical  image  reconstruction  for  biomedical  imag¬ 
ing  and  diagnostics  has  attracted  considerable  attention  due  to  in  part  its  potential  for 
noninvasive  clinical  applications.1-3  Near-infrared  light  can  probe  the  internal  struc¬ 
ture  (absorption  and  scattering  coefficients)  of  a  highly-scattering  turbid  medium  such 
as  human  tissues  by  measuring  scattered  light  around  the  medium.  The  propagation 
of  light  in  turbid  media  is  governed  by  the  radiative  transfer  equation.4  Approxima¬ 
tions  based  on  truncation  in  the  spherical  harmonic  expansion5  or  cumulant  expansion 
of  the  photon  distribution  function6  7  are  used  in  image  reconstruction.  However,  the 
numerical  reconstruction  is  extremely  complicated  and  computational  expensive  even 
after  approximation  and  linearization  of  the  forward  problem. 

One  of  the  recent  developments  in  optical  tomography  is  to  apply  Fourier  tech¬ 
nique  to  analyze  the  light  propagation  in  turbid  media  (so  called  “optical  diffraction 
tomography”).  Diffraction  tomography  may  be  traced  back  to  Emil  Wolf’s  original 
work  of  determination  of  the  spatial  distribution  of  the  refractive  index  from  measure¬ 
ments  of  the  intensity  transmission  functions  of  holograms.8’ 9  Diffraction  tomography 
has  been  extensively  applied  to  imaging  using  ultrasound  and  microwave.10-12  Re¬ 
cently  there  are  considerable  effort  to  apply  this  diffraction  technique  to  tomographic 
imaging  using  diffuse  photons;13-15  in  particular,  Schotland  and  Markel’s  work.1617 
The  benefit  of  the  diffraction  formalism  of  optical  imaging  comes  from  the  diago- 
nalization  of  the  point  spread  function  (the  Green’s  function)  in  the  spatial  Fourier 
space  in  presence  of  a  translational  invariance,  yielding  a  significant  simplification  in 
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reconstruction. 


Since  light  is  highly  scattered  inside  the  turbid  medium  and  very  different  spatial 
distributions  of  the  optical  property  in  the  medium  may  result  in  indistinguishable 
observations,  the  reconstruction  of  the  internal  structure  from  scattered  light  based  on 
a  forward  model  formulated  in  either  the  real  space  or  the  Fourier  space  (diffraction 
tomography)  is  not  unique  and  hence  ill-posed.18  The  inevitable  photon  counting 
noise  accentuates  this  problem  further.  A  common  practice  is  to  add  regularization 
to  stabilize  the  inversion  process. 

One  of  the  most  popular  regularization  methods  was  proposed  by  Tikhonov.19 
Describe  an  imaging  of  a  set  of  measurements  b  due  to  an  internal  structure  f  of  a 
medium  as: 

At  4-  n  =  b  (1) 

in  a  framework  of  a  linearized  reconstruction  where  n  is  the  noise  involved  in  measure¬ 
ment.  Tikhonov  seeks  a  best  possible  solution  to  Eq.  (1)  by  solving  an  optimization 
problem 

min{||Af-b||2  +  A2||Lf||2}  (2) 

which  balances  the  goodness  of  fit  (the  first  term)  and  the  closeness  to  the  prior  infor¬ 
mation  (the  second  term).  Here,  A  is  the  regularization  parameter  which  determines 
the  trade-off  and  L  is  a  penalty  operator  which  draws  the  solution  toward  its  null 
space  (chosen  to  approximate  the  prior  information).20 

An  assumption  of  white  noise  in  all  the  observations  at  one  common  level  is 
implicit  in  Eq.  (2)  whereas  the  noise  presented  in  photon  counting  tends  to  be  Poisson 
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and  position-dependent,  proportional  to  the  square  root  of  light  intensity  at  each 
detecting  element  in  the  shot- noise  limit.  Intuitively,  regularization  should  consider 
both  the  noise  presented  in  observation  and  the  prior  information  of  the  object. 

We  will  show  in  this  paper  a  proper  modeling  of  the  noise  improves  the  quality 
of  reconstruction.  A  generalized  Tikhonov  regularization  formalism  is  used  in  image 
reconstruction  to  incorporate  both  the  prior  information  and  the  noise.  Furthermore, 
the  optimal  regularization  parameter  will  be  shown  linked  directly  to  the  prior  signal 
to  noise  ratio  (PSNR). 

This  paper  is  organized  as  following.  We  will  first  review  the  diffraction  tomog¬ 
raphy  using  diffuse  photons.  We  present  an  explicit  regularized  inversion  formula  for 
diffraction  image  reconstruction  after  discussion  of  a  generalized  Tikhonov  regular¬ 
ization  method.  Computer  simulation  results  are  shown  to  demonstrate  improvement 
of  the  quality  of  reconstruction  by  the  new  inversion  formula  using  time-resolved 
data  with  realistic  noise  added  for  reconstruction  of  a  planar  geometry  using  multiple 
sources  and  detectors. 

2.  Theory 

A.  Optical  diffuse  diffraction  tomography 

The  formalism  for  diffuse  diffraction  tomography  has  been  given  by  various  authors 
in  time-  and  frequency-  domains.13  15  17  Follow  the  Dirac  notation  used  by  Markel 
and  Schotland,16  diffraction  formalism  of  the  optical  imaging  using  diffuse  photons  is 
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reviewed  here.  The  diffusion  equation  in  time  domain  is 


dt\u(t))  +  H\u(t))  =  \S(t))  (3) 

where  |u)  is  the  energy  density,  the  Hamiltonian  H  =  -V  •  D(r)V  +  a(r),  D(r) 
and  a(r)  are  the  position-dependent  diffusion  and  absorption  coefficients,  and  |S)  is 
the  source  power  density.  This  is  the  Schrodinger  equation  in  imaginary  time  where 
ct(r)  can  be  interpreted  as  the  interaction  potential  and  D(r)  as  the  inverse  of  the 
position-dependent  mass.  The  evolution  operator  (the  Green’s  function)  is  given  by 

G  =  0(f)  exp(-Ht)  (4) 

where  0(f)  is  the  Heaviside  step  function. 

It  is  best  to  regard  the  diffusion  equation  (3)  as  an  equation  in  a  Hilbert  space  and 
a  basis  of  position  r  and  time  t  is  used  in  Eq.  (3).  Some  other  commonly  used  basis 
include  the  temporal  frequency  oj  and  the  lateral  spatial  frequency  q.  The  rotation 
from  one  basis  to  another  can  be  facilitated  using  the  following  identities 

J\t)dt(t\  =  J\u)^(u>\  =  J\p)(Pp(p\  =  j  |q)AS_(q|  =  l  (5) 

(w|t)  =  exp (iujl),  (q|p)  =  exp(-iq  ■  p) 

due  to  the  closure  property  of  the  respective  basis. 

In  one  class  of  problems,  the  position-dependent  diffusion  and  absorption  coeffi- 
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dents  D( r)  and  a(r)  can  be  decomposed  to  D( r)  =  D0+<5D(r)  and  a(r)  =  ft0+<Sa(r) 
where  5D( r)  and  Sa(t)  can  be  regarded  as  a  perturbation  to  the  values  D0  and  a0  of  a 
uniform  background.  The  scattered  field,  i.e.,  the  perturbation  of  the  energy  density, 
due  to  the  embedded  inhomogeneities  5Z)(r)  and  <5o:(r),  is  given  by 


<t>s  =  <j>-<j>o  =  G\u)-G0\u)  =  -G0AHGo\u)  (6) 

to  the  first  order  Born  approximation  where  ©o  is  the  energy  density  from  unperturbed 
uniform  background,  6  is  the  total  energy  density  from  the  perturbed  medium,  and 
Aff  =  -V-&D(r)V  +  <5a(r). 

Let’s  restrict  ourselves  to  a  planar  geometry  which  possesses  a  translational  invari¬ 
ance  along  the  lateral  directions.  The  Green’s  function  of  a  uniform  planar  geometry 
has  a  form  of 


<p2(|G„|pV('>  =  &(t-f)gs-L-—ac p 


r  (p  -  p')2 

4 D0(t  - 1') 


-  a0{t  - 


Gz(z,z',t-t') 


(7) 


where  Gz  is  the  axial  component  of  the  Green’s  function  and  dependent  on  the  specific 
boundary  condition.  Its  diffraction  counterpart  is  given  by 


(qzt\G o\q  z't')  =  J  (q|p)  d2 p  (pz't\G0\p'z't')  d2p'  (p'\q')  (8 

=  e(^  “  0  exp  \-{Dq 2  +  a){t  -  t') }  Gz(z ,  z',  t  -  t'){2n)2S(q  -  q') 
=  G0(q.z.z,t  -  t')(2Tr)25(q  —  q') 
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by  changing  the  basis  p  to  q  in  time  domain.  The  benefit  of  a  diffraction  formalism 
with  the  lateral  spatial  frequency  q  is  now  clear  that  the  Green’s  function  is  diagonal 
with  respect  to  the  basis  and  has  a  most  simple  formulation. 

The  scattered  field  is  then  expressed  as: 


=  -jf  dT  J  dzGo(qd,zd,z,t-T){qdZT\AH\qsZT)G0(qS:z,zs,T-t0 ) 

(9) 


from  Eq.  (6)  using  the  identities  (5)  due  to  an  incident  “point”  source  |q.,2s)  6(t-t0). 
The  operator  {qdZT\AH\qszr)  can  be  further  evaluated  in  a  similar  fashion  using  (5). 
This  results  in15 


q.s,  t  to)  —  J  dz  [^(q^,  q5j  z)5a{q<i  —  qfi,  z)  +  ^(q^,  q5?  z)8D(qcI  —  q5,  z)] 

(10) 

where  the  weight  function  for  absorption  and  diffusion  inhomogeneities  is  given  by 


ka 

kd 


I  8/TG(c^fji t  Zdi  Z)  t  T)G’(qS7  Z,  T  to) 
Jto 


f 


dr 


f dGjcy,  z_g  z  t  -  t)  dGjgs,  z^t  -  t0) 


to  dz  dz 

+qs  •  qdG(qd,  zd,  z,t  -  r)G(qs,  z,  zs,t  -  t0)] 


(11) 


After  performing  a  variable  change, 


q<f  =  P  +  q/2 
qs  =  p-q/2 
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the  forward  problem  can  be  rewritten  as  an  integral  equation 

<k(P  +  q/2,  P  -  q/2 )  =  j  dz  [ka( p,  q,  z)5a(q,  z)  +  /c^(p,  q,  z)5D{ q,  z)}  (13) 

up  to  a  constant  multiplier  (—1  here).  The  equation  (13)  defined  an  independent 
integral  equation  of  5o:(q,  z)  and  5D( q,  z)  for  each  fixed  q  in  the  form  of  Eq.  (1).  The 
corresponding  result  in  frequency  domain  is  obtained  replacing  the  integral  over  time 
by  a  factor  2ir5(u!d  —  u.'s)  and  changing  t  —  r  and  t  —  to  to  and  u>$  respectively, 
which  had  been  given  in  Markel  and  Schotland.16 

The  formal  inversion  formula  for  Eq.  (13)  with  a  fixed  q  is  given  by 

6a(q,z)  =  J  d2p'dz'  (z|T~1(q)|2/)  K*A(p' ,  q,  z')<p(pr  +  q/2,  p'  —  q/2)  (14) 

SD(q,  z)  =  J  d2p'dz’  (z\T~\q)\z')  K*D(p\  q,  z')4>{ p'  +  q/2,  p'  -  q/2) 

where 

(z|T(q)|z')  =  J d2p [k^(p, q,  z)/c/1(p, q,  z')  +  K^(p, q,  z)kd(p, q,  z')\ .  (15) 

This  can  be  verified  by  substitution  of  Eq.  (14)  into  Eq.  (13).  This  is  the  Moore- 
Penrose  generalized  inverse  of  the  linear  system  (13).21 

The  incident  Fourier  ‘point”  source  |q szs)  6{t  -  ts)  in  time  domain  or  Iq,-;,^)  in 
frequency  domain  can  be  generated  from  a  spatially  modulated  mask  under  illumi¬ 
nation  of  a  plane  wave,  or  by  Fourier  transform  of  an  array  of  point  sources. 
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A  scheme  is  given  in  Fig.  (1)  where  a  point  source  scans  grid  by  grid  covering 
the  input  window  and  the  emitted  light  on  the  exit  window  is  captured  by  a  CCD 
camera  in  a  planar  geometry. 

We  proceed  to  provide  a  regularized  version  of  the  above  solution  (14)  after 
introduction  of  a  generalized  Tikhonov  regularization  method. 

B.  Generalized  Tikhonov  regularization 

The  inverse  problem  of  the  form  of  Eq.  (1)  is  ubiquitous  and  regularization  is  usually 
required  to  stabilize  the  inverse  process.  As  we  have  pointed  out  in  the  introduction, 
the  standard  Tikhonov  regularization  scheme  assumes  implicitly  an  uniform  distribu¬ 
tion  of  white  noise  in  all  observations  which  contradicts  with  the  position-dependent 
noise  of  photon  counting  from  detecting  elements  in  optical  imaging. 

An  intuitive  way  of  dealing  the  inverse  problem  Eq.  (1)  is  to  regard  f.  n  and 
b  as  a  realization  of  Gaussian  weak  random  variables  /,  n  and  b  in  the  Hilbert 
space  L2(R).22  A  Gaussian  weak  random  variable  is  uniquely  determined  by  its  mean 
element  and  covariance  operator.  Under  the  assumptions  that  f  and  n  have  zero  mean 
and  uncon elated,  and  that  the  covariance  operator  Rnn  has  a  bounded  inverse,  the 
best  linear  estimate  of  / ,  namely  Cob,  which  minimizes  the  mean  square  error,  was 
shown  by  Viano  et.  al.23  as 


C0  =  RffA*(ARffA*  +  Rnn)~'  (16) 

where  a  superscript  "*"  means  hermitian  and  the  covariance  is  defined  as  Rr„  = 

.cy 
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E  {[x  —  E(x)]  \y  —  E(y)]  }.  This  formula  is  derived  by  minimizing  the  mean  square 
error  and  closely  connected  with  the  Kalman  filter  theory.24  It  is  easy  to  verify  that 
Co  can  be  rewritten  in  a  more  familiar  form: 

Co  =  (A*R£A  +  (17) 

Let’s  introduce: 

'  Rnr,  =  crV.  trace(V')  =  1  (IS) 

RfJ  =  t2(L*L)~\  trace((L*T)-1)  =  1 

where  <r  and  t  can  be  interpreted  as  the  magnitudes  of  the  noise  and  the  prior  signal 
respectively.  Eq.  (17)  yields 

C0  =  (A*V~lA  +  ^L*L)~1A*V~1.  (19) 

The  linear  estimate  Cob  is  now  recognized  as  the  solution  of  a  generalized  Tikhonov 
regularization  problem:25 

cugmin  ^ (Af  -  b)*y-J(Af  -  b)  +  ^  |jLf||2'*  .  (20) 

Moreover,  the  choice  of  the  regularization  parameter  A  =  cr/r  must  be  optimal  for 
this  inverse  problem  with  the  prescribed  noise  covariance  and  prior  information  (18) 
since  Cob  is  the  best  estimator.  The  standard  Tikhonov  regularization  responds  to 
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the  case  where  V  —  I  where  I  is  an  identity  matrix.  The  penalty  term  is  commonly 
set  also  to  an  identity  matrix  in  the  standard  Tikhonov  regularization.^ 

This  shows  that  the  inverse  of  the  optimal  regularization  parameter  A*: 

1/A*  =  PSNR  =  rja  = 

is  linked  directly  to  the  prior  signal  to  noise  ratio  (PSNR).  The  value  of  the  optimal 
regularization  parameter  is  proportional  to  the  magnitude  of  the  noise  and  inverse 
proportional  to  the  magnitude  of  the  prior  signal.  Note  PSNR  is  different  from  the 
conventional  signal  to  noise  ratio  (SNR)  which  usually  means  the  ratio  of  magnitudes 
of  noiseless  observation  to  noise. 

C.  Noise  analysis  and  regularized  diffraction  reconstruction 

Measurements  captured  by  the  CCD  camera  on  the  exit  window  is  noise  contami¬ 
nated.  The  noise  is  Poisson  with  its  variance  given  by  Var  [n(pd,  ps)]  =  b(pd,  ps)  in 
the  shot  noise  limit.  Here  we  have  assumed  the  detector  element  has  a  unit  size  one 
and  b(pd,  pA)  is  the  count  of  photons  on  the  pd  detecting  element  from  a  point  source 
at  ps .  The  covariance  of  the  noise  in  our  Fourier  space  is  thus 

-^ri(qrf.qs)n(q^q')  =  E  ^(q(/|n|qs)  (q^|n jq^)j  -  E  [(qdHqs)]  E  (q^|n|q^)j  (22) 

=  J  d2pdd~ps  exp  [-i( qd  -  q Q-  pd  +  i(q$  -  q()  ■  ps]  Var  [n(pd.  p5)| 
=  (qd  -  q'iMq*  -  q's) , 


trace  (R//) 
trace  ( Rnn ) 
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assuming  that  the  noise  at  different  pixels  and  from  different  sources  is  independent. 

One  important  consequence  of  this  noise  covariance  (22)  is  that  the  noise  covari¬ 
ance 

#nn(q)  =  #n(q,.qs)n(q^),  where  qd  -  qs  =  -  q's  =  q  (23) 

in  one  subspace  of  a  fixed  q  is  circulant  and  symmetric.  This  yields 

(pl-Rnn(q)lp')  =  j  (p|p)  d2p  (p|A|p')  d2p'  (p'  |p')  (24) 

where 

;  # 

(p|A|pr)  =  Sb(p,  0)5(p  —  p')  (25) 

is  diagonal.  The  inverse  of  i?nn(q)  is  then  simply  given  by 

(pK(q)lp')  =  J  (p|p>  d2p  (p\l\-l\p')  d2 p'  (pV)  (26) 

where 

(PIA’V')  =  S~lb-\p,Q)5{p  -  p').  (27) 

This  can  be  easily  verified  by  evaluating  /  <p|HBW(q)|p'>  0  <p'|iCKq)|p")  =  <p|p">- 
The  generalized  Tikhonov  reconstruction  of  our  formal  solution  (14)  can  now  be 


T2 


regularized  to 


<fo(q,  z) 
5D(q,  z ) 


/  (*I^S(q)k)  «a(p',  q,  *')  (p'liCn (q)|p">  4>{ p"  +  q/2,  p"  -  q/2) 

/ d2p'd2p"dzf  (z\T-}(q)\z')  K*D(p ',q,z')  (p'|^(q)|p")  <£(p"  +  q/2,  p"  -  q/2) 

(28) 


from  Eq.  (19)  where 

{2|Treg(q)|,::')  =  J  d2pd2p'  \ka(p, q,  z)  (p|R~n(q)|p')  ka(p',  q,  z')+  (29) 

«r>(P,  q>  z)  (p|i?^(q)|pO  K£>(p',q,z')l  +  (zl^J^q)^') 

and  (p|-R^(q)|p')  is  given  by  Eq.  (26). 

Image  reconstruction  is  carried  out  by  first  obtaining  <5a(q,2)  and  5D(q,z)  for 
each  q  by  Eq.  (28)  formulated  in  subspace  q.  The  perturbations  in  absorption  and 
diffusion  coefficients  in  real  space  are  obtained  from  an  inverse  Fourier  transform 
afterwards. 

D.  Prior  information  and  regularization 

The  prior  signal  covariance  Rff,  or  equivalently,  the  penalty  term  L  in  Eq.  (18),  should 
be  chosen  according  to  the  prior  information.  For  simplicity,  we  will  only  consider  the 
absorption  inhomogeneity  5n( r)  here.  If  the  prior  signal  <5a(r)  is  assumed  to  be  not 
correlated  between  different  voxels  in  real  space,  corresponding  to  a  choice  of  L  =  I 
where  I  is  an  identity  matrix,  the  prior  signal  covariance  in  the  Fourier  space  is  then 
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given  by: 


—  5{z  —  z')  J  d2pexp  [— i(q  —  q')  •  p)\  Var  [5oc(p,  z)\ , 


analog  to  the  derivation  of  the  noise  covariance  (22).  The  prior  signal  to  noise  ratio 
in  a  subspace  q  is  found  to  be 


PSNR  =  5-.  //jtVarMr)l 

V  I  d2pb{p,  o) 


(31) 


independent  of  q  where  S  gives  the  dimensionality  (the  total  number  of  values  of  p) 
of  the  subspace. 

This  mandates  a  common  regularization  parameter  to  be  used  in  all  subspaces  q 
from  the  optimal  regularization  criterion  (21).  As  the  magnitude  of  the  prior  signal 
is  usually  unknown,  additional  regularization  parameter  choice  rules,  such  as  the  L- 
curve  method,26  must  also  be  invoked.  The  regularization  parameter  in  this  paper  is 
obtained  by  locating  the  L-corner  of  the  global  L-curve  (£q  ||/?(q)||2  vs  £q||f?(q)||2 

where  p(q)  and  r}(q)  are  the  residue  of  error  and  the  size  of  solution  respectively  in 
the  subspace  q). 

The  generalized  Tikhonov  regularization  does  not  restrict  the  prior  signal  covari¬ 
ance  (or  the  penalty  term)  to  be  an  identity  matrix.  Some  other  forms  of  the  penalty 
term  is  also  commonly  used,  including  a  matrix  of  suitably  weighted  first  or  second 
order  differences  for  a  prior  smooth  signal25  and  spatial  variant  version.27  In  this  pa¬ 
per,  we  will  test  both  the  identity  and  the  first  order  differential  penalty  in  image 


reconstruction. 


3.  Results 

To  demonstrate  our  regularized  solution,  we  perform  computer  simulations  using  a 
time-resolved  scanning  source  scheme  [see  Fig.  (1)].  Without  losing  generality,  we  will 
assume  a  uniform  distribution  of  diffusion  coefficient,  i.e. ,  8D  —  0.  In  our  simulations, 
the  input  and  exit  windows  are  divided  into  21  x  21  grids  and  20  grids  along  the  axial 
direction.  Each  voxel  has  a  volume  of  3  x  3  x  1.5mm^.  The  distance  between  the 
input  and  exit  windows  of  the  slab  is  30mm.  The  medium  assumes  a  refractive  index 
1.38,  an  absorption  coefficient  0.003mm_1  and  transport  mean  free  path  1.25mm, 
comparable  to  those  found  in  a  typical  human  breast. 

Two  absorptive  inhomogeneities  are  located  at  (14, 14, 10)  voxel  with  100%  in¬ 
crease  in  absorption  than  that  of  background  and  at  (5, 5, 15)  voxel  with  50%  increase 
in  absorption  than  that  of  background.  Poisson  noise  is  added  to  the  simulated  photon 
counts.  Note,  the  added  noise  is  proportional  to  the  square  root  of  the  total  photon 
counts  (p o  +  <t>s  with  the  absorbing  objects  embedded  and  not  the  scattered  field  4>s. 
Only  snapshots  of  photon  counts  at  time  delay  of  400ps  are  used  in  the  image  recon¬ 
struction  here  as  we  are  most  interested  to  see  the  effect  of  modeling  of  noise  on  the 
quality  of  image  reconstruction.  Inclusion  of  different  time  delays  is  beneficial  to  the 
reconstruction. 

To  quantify  the  noise  added  in  the  simulation,  we  define  the  signal  to  noise  ratio 
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(SNR)  and  the  effective  signal  to  noise  ratio  (SNR’)  in  our  simulation  as 


SNR(dB)  =  10  log10 
SNR'(dB)  =  10  log10 


2Zps  11*^0 (pdj  Ps)  &s(,Pdi  Ps)|| 

EPdEPs  ll«(pd,Ps)ll2 
EPdEps  IIMp^p*)!!2 
EprfEPs  \Hpd,ps)\\2 


(32) 


where  (f>o(pd ,  ps)  is  the  noiseless  observation  at  Pd  from  a  point  source  at  ps  in  the 
uniform  background  and  <j>s{pd,  Ps)  is  the  perturbation  of  this  observation  due  to  the 
existence  of  the  inhomogeneities.  Our  simulation  is  performed  with  SNR  =  80dB, 
and  an  effective  SNR'  =  4.4dB.  This  noise  is  significant  as  the  total  amount  of  the 
scattered  field  is  less  than  that  of  noise.  One  sample  profile  is  given  in  Fig.  (2)  to 
show  the  degree  to  which  the  scattered  field  <z>,  is  affected  by  the  added  Poisson  noise. 

The  reconstruction  result  using  the  standard  Tikhonov  regularization  is  shown  in 
the  following  figures.  Fig.  (3)  gives  the  L-curve  (the  residue  of  error  ||/>||2  =  \\Ax  -  b||2 
vs  the  size  of  solution  ||t?||2  =  ||x||2).  The  L  corner  suggests  an  optimal  regularization 
parameter  to  be  about  0.2.  This  turns  out  to  agree  with  the  inverse  of  the  prior  signal 
to  noise  ratio  PSNR  =  0.2  calculated  using  Eq.  (31).  The  reconstruction  result 
using  this  regularization  parameter  A  =  0.2  is  given  in  Fig.  (4).  The  two  inhomo¬ 
geneities  at  (14, 14, 10)  and  (5, 5, 15)  are  resolved  ■with  a  correct  central  position  and 
minimal  blurring  over  the  transversal  direction.  Their  axial  position,  however,  is  not 
resolved  as  well  as  in  the  transversal  direction.  The  axial  central  position  of  the  first 
inhomogeneity  (14, 14, 10)  is  resolved  to  be  on  the  right  layer  while  the  second  one 
(5, 5, 15)  peaks  on  the  layer  17  rather  than  the  right  one  (layer  15).  Both  resolved  in- 
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homogeneities  expand  considerably  over  the  axial  direction.  As  a  result,  the  resolved 
absorption  coefficients  are  only  a  fraction  of  the  input  values  (100%  and  50%  are 
expected  for  the  first  and  second  inhomogeneities  respectively). 

The  reconstruction  result  for  the  same  data  using  the  generalized  Tikhonov  reg¬ 
ularization  with  an  identity  penalty  (L  =  I)  is  shown  in  the  following  figures.  Fig.  (5) 
gives  the  L-curve  (the  residue  of  error  ||p||2  =  (Ax  -  b)*F'1(Ax  -  b)  vs  the  size  of 
solution  ||t?||2  =  ||x||2).  The  optimal  regularization  parameter  suggested  by  the  L  cor¬ 
ner  agrees  with  PSNR-1  again,  both  pointing  to  A  =  0.2.  The  reconstruction  result 
using  this  optimal  regularization  parameter  is  given  in  Fig.  (6).  Comparing  to  the 
reconstruction  given  by  the  standard  Tikhonov  regularization  in  Fig.  (4),  noticeable 
improvement  is  observed  in  the  resolution  of  the  axial  position  of  the  second  inhomo¬ 
geneity  which  peaks  now  around  layers  16  and  17.  closer  to  the  right  value  (layer  15). 
The  ratio  between  the  resolved  absorption  coefficients  of  the  two  inhomogeneities  is 
also  closer  to  that  of  the  input  values. 

A  second  generalized  Tikhonov  reconstruction  is  performed  with  a  first  order 
differential  operator  penalty  on  the  same  data.  Fig.  (7)  gives  the  L-curve  (the  residue 
of  error  ||p||2  =  (Ax  -  b)*l/-1(Ax  -  b)  vs  the  size  of  solution  |j7?||2  =  ||Lx||2).  This 
L-curve  is  not  as  well  formed  as  the  previous  two  other  L-curves  given  in  Figs.  (3) 
and  (5).  The  regularization  parameter  A  =  0.2  is  shown  on  the  L-curve  figure  which 
is  now  to  the  left  of  the  L-corner.  The  reconstruction  result  using  the  regularization 
parameter  A  =  0.2  is  given  in  Fig.  (8).  The  axial  position  of  the  two  inhomogeneities  is 
now  correctly  resolved.  Narrower  axial  profiles  are  observed  for  both  inhomogeneities 
(widths  at  half  peak  of  8  layers  for  both  inhomogeneities  compared  to  a  width  of 
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half  peak  of  12  layers  for  the  first  inhomogeneity  and  a  half  width  of  peak  of  7 
layers  for  the  second  inhomogeneity  using  the  generalized  Tikhonov  regularization 
with  an  identity  penalty).  In  addition,  the  resolved  absorption  coefficients  are  more 
closer  to  the  input  values  (peak  values  of  0.115  and  0.108  are  obtained  for  the  two 
inhomogeneities  respectively). 

A  better  reconstruction  can  be  obtained  by  using  an  aggressive  regularization 
parameter  A  =  0.02  at  the  expense  of  a  degraded  smoothness  of  the  reconstruction 
[see  Fig.  (9)].  In  addition  to  the  correct  resolution  of  the  lateral  and  axial  position  of 
both  inhomogeneities,  the  resolved  absorption  coefficients  over  that  of  the  background 
boc/a o  now  approach  0.31  and  0.19  for  the  first  and  second  inhomogeneities  respec¬ 
tively,  the  best  among  all  the  reconstructions.  These  resolved  peak  values  correspond 
to  31%  and  38%  of  the  input  values  respectively. 

4.  Discussion 

Our  results  of  computer  simulation  of  reconstruction  inside  a  planar  geometry  demon¬ 
strated  the  importance  of  modeling  correctly  the  noise  presented  in  measurements. 
With  a  proper  modeling  of  the  Poisson  noise  and  an  appropriate  regularization,  a  bet¬ 
ter  axial  localization  and  optical  property  reconstruction  is  obtained  in  our  simulation 
with  realistic  noise  added. 

The  performance  of  the  standard  Tikhonov  regularization  is  not  as  good  as  the 
generalized  version  but  not  so  bad  in  the  diffraction  tomography  thanks  to  the  special 
form  of  the  noise  covariance  Rnn( q)  in  a  subspace  q.  The  most  significant  elements  of 
the  noise  covariance  Rnn{ q)  given  in  Eq.  (24)  is  along  the  diagonal  and  of  the  same 
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value.  The  Fourier  transform  of  the  observation  data  in  a  diffraction  tomography  in 
effect  whitens  the  noise.  This  hence  allows  an  approximation  of  the  noise  covariance 
in  the  form  of  an  identity  matrix,  which  happens  to  be  used  by  the  standard  Tikhonov 
regularization.  Optical  imaging  using  diffuse  photons,  in  general,  needs  to  consider 
the  noise  explicitly. 

The  appropriate  regularization  for  an  inverse  problem  is  an  art  rather  than  sci¬ 
ence.  The  linkage  between  the  optimal  regularization  parameter  and  the  prior  signal 
to  noise  ratio  (PSNR)  may  provide  valuable  guidance  in  the  choice  of  the  regulariza¬ 
tion  parameters.  The  fact  that  the  optimal  regularization  parameter  is  proportional 
to  the  noise  magnitude  and  inverse  proportional  to  the  prior  signal  is  consistent  with 
intuitions.  Even  in  the  absence  of  known  prior  signal  which  is  usually  the  case,  this 
linkage  helps  to  choose  the  optimal  regularization  parameter.  In  our  computer  simu¬ 
lation,  all  the  regularization  parameters  in  different  subspace  q  are  set  to  the  same 
value  according  to  this  criterion.  We  have  also  observed  that  a  first  order  differen¬ 
tial  penalty  performs  better  than  an  identity  penalty  in  the  generalized  Tikhonov 
regularization  method  for  optical  image  reconstruction  using  diffuse  photons. 

The  explicit  inversion  formula  developed  in  this  paper  is  different  from  that  pro¬ 
posed  by  Markel  and  Schotland.16  We  formed  our  inversion  formula  based  on  a  gener¬ 
alized  Tikhonov  regularization  formula  (19).  We  reduced  the  image  reconstruction  to 
a  series  of  one  dimensional  inverse  problem  of  the  dimension  of  the  number  of  discrete 
divisions  in  the  axial  direction,  compared  to  a  series  of  two  dimensional  inversion  of 
a  plane  with  a  fixed  axial  location  used  in  Ref.  16.  The  number  of  operations  re¬ 
quired  inci  eases  only  linearly  with  the  total  number  of  detector-source  pairs  in  our 
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formulation  and  can  be  easily  extended  to  use  the  cumulant  approximation  of  photon 
migration  in  turbid  media.6’7  Another  strength  is  the  flexibility  in  choosing  a  suitable 
penalty  term.  We  also  want  to  point  out  that  the  inversion  formula  in  Ref.  16  can 
be  enhanced  to  include  the  noise  and  prior  information  terms  using  the  best  linear 
estimator  given  in  Eq.  (16),  which  will  be  studied  later. 

In  summary,  an  explicit  regularized  inversion  formula  for  a  three-dimensional 
diffraction  tomographic  imaging  using  diffuse  photons  is  derived.  This  approach  in¬ 
corporates  both  the  prior  information  and  the  noise  explicitly  using  a  generalized 
Tikhonov  regularization  scheme  in  which  the  optimal  regularization  parameter  links 
directly  to  the  prior  signal  to  noise  ratio.  Proper  modeling  of  noise  and  appropriate 
regularization  is  shown  to  improve  the  quality  of  image  reconstruction  of  the  optical 
diffraction  tomography. 
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